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ABSTRACT 

A binary collision simulation of the radiation damage problem is 
proposed. Energy thresholds and parameter adjustment are discussed. 
Preliminary results show strong focusing in the (110) direction and 
weaker focusing in the (100) direction. 

Although not complete, the model is expected, with suggested 
modifications, to adequately represent the radiation damage problem. 
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of the U. S. Naval Postgraduate School in this investigation. The authors 
also are pleased to acknowledge the assistance of Lt. Herbert L. Hoppe, 
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the development of the binary collision interaction subroutine. 
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1. Introduction 
The interaction of energetic particles with matter may be divided 
into two energetic regions of interest: 

1) High energy events in which interactions are primarily 
inelastic electronic processes of excitation and ionization, 
and 

2) Low energy events in which the interactions are essentially 
elastic interactions between whole atoms. 

There is an intermediate energy range in which both processes are im- 
portant, but until the two limiting cases are better understood, attempts 
to understand the intermediate range seem futile. i* The threshold energy 
for which electronic excitation and ionization processes become impor- 
tant depends upon the mass of the energetic particle. An approximate 
formula for this threshold energy is E=A kev., whereA is the atomic 
mass of the energetic particle. The threshold for electronic processes 
seems to be relatively insensitive to the atomic mass of the target | 
acer For metals this threshold is of the order of 10-100 kev. The 
high energy processes of electronic excitation and ionization have re- 
ceived considerable attention and are reasonably well understood. 3-6 
We will confine our attention to the energy region below this threshold 
for inelastic electronic processes. 

Many physical properties of a material are changed when the 
material is aempiea. The Young's modulus and the electrical 
* All footnotes refer to the BibNogeEpty. 
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resistivity are among those properties of materials which show most 
significant changes. For metals and semiconductors, the residual 
electrical resistivity is most sensitive to radiation. 

Semiconductors, for example, may show changes in electrical 
resistivity of up to four orders of magnitude if but there is no known 
simple relation between the magnitude of the change in electrical 
resistivity and the amount of radiative energy absorbed by a semi- 
conductor. Radiation changes some semiconductors from n-type to the 
p-type. 

The radiation damage problem involves chemical bonding between 
the atoms of the target material as well as elastic collisions between 
atoms. For low energy collisions, chemical bonding forces vary widely 
for different atomic species, but are of the same order of magnitude as 
the elastic repulsive forces. Therefore, to reduce the complexity of 

the problem, we have limited our investigation to simple crystals of 
face centered cubic metals, with special emphasis on copper. 

Davies, et al. have developed a chemical technique for removing 
thin layers (= 37 A) from aluminum gouteee The thickness of the layers 
is easily controlled. The technique has been widely used to study the 
range of radioactive ions in aluminum.9 «19 Ranges have been measured 
for Csl37, Rb&6 | Na24, and K42, For example, the mean range of 
10.5 kev. Na24 ions in aluminum was found to be about 250 Angstroms. 

When a metal is irradiated with heavy energetic ions (>_50 kev.) a 
significant number of atoms are ejected from the target through the surface 
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which is bombarded by the incident ion beam. This phenomenon is 


known as sputtering. Under suitable conditions as many as 40 atoms or 


ll 


ions may be sputtered for each incident particle. Examination of the 


angular distribution of sputtered particles reveals definite preferential 
ejection in directions which correspond to (110), (100), and (111) lattice 
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2. Lattice Effects 

Whether a lattice atom leaves its original site or not after a 
collision depends upon the energy absorbed by the atom and its final 
direction of motion with respect to the lattice. If an energetic lattice 
atom does not leave its site, it will transfer energy to its neighbors 
through elastic collisions until thermal equilibrium is reached. Be- 
cause of crystalline structure, the dissipation of energy is expected to 
be anisotropic. The only exception to the eventual thermalization of 
the excess energy of a bound lattice atom is the case in which sufficient 
energy reaches a surface of the material to knock one or more atoms away 
from the surface. 

If a lattice atom is given sufficient energy to leave its site per- 
manently, and no other atom replaces the energetic atom, a vacant 
lattice site is created. The Surrounding lattice atoms relax slightly 
toward the vacant site, but the lattice structure remains essentially the 
same as that of a perfect lattice. 17 The energetic atom travels through 
the lattice giving up energy through elastic collisions with the atoms 
near its path until it no longer has sufficient energy to force its way 
between the atoms of the crystal. If there are no vacant lattice sites 
near the position at which an energetic atom is brought to rest, it will 
be forced to occupy some position other than a regular lattice site. 

Three possible equilibrium positions have been suggested for these 
interstitial atoms in a face centered cubic crystal: the crowdion, the 
split interstitial pair, and the center of a unit cell. 7 The crowdion 
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consists of an extra atom in a (110) line of atoms. In the split inter- 
stitial pair configuration, the interstitial shares a regular lattice site 
with another lattice atom so that each is symmetrically displaced from 
the equilibrium position in a (100) direction. An atom in the center of 
a unit cell may also be described as a crowdion in a (100) direction. 
The lattice atoms surrounding an interstitial are displaced from their 
regular lattice sites to accommodate the extra atom. More will be said 
about the interstitial configuration later. 

An energetic atom (knock-on) may transfer essentially all of its 


18 In this case, 


energy to a single lattice atom in a head-on collision. 
the lattice atom becomes the knock-on, leaving a site for the former 
knock-on to occupy. The new knock-on may, in turn, transfer essentially 
all of its energy to another atom. Thus, a chain of replacements may 
transport an interstitial to a position which is several times the dis- 
tance between nearest neighbors from the corresponding vacancy. Re- 
placement chains may be cyclic so that the last replacement in the chain 
fills the initial vacancy, leaving no vacancies or interstitials. ms The 
configuration consisting of a vacancy and a nearby interstitial, known 

as a Frenkel pair, may be stable or unstable. The stability of a Frenkel 
pair depends upon both the distance of separation and the orientation 
(with respect to the lattice) of the interstitial and the vacancy. 


Silsbee has Showa? 


that an isolated line of hard spheres may 
focus and transfer energy along the line. The condition for focusing 


is that s/d< cos 6 , where s is the distance between adjacent spheres, 


3 





d is the diameter of the spheres, and 6 is the angle between the 
velocity of the first atom in the line and the line. The hard core 
approximation (diameter of spheres equal to distance of closest 
approach for a head-on collision) is useful in this work. For reason- 
able repulsive potentials such as the Coulomb, Bohr, Born-Mayer, etc., 
this modified model predicts a focusing effect below a threshold energy, 
and a de-focusing effect above this threshold energy. 17 When the 
Gibson potential #2 simulates copper, this model predicts focusing in 
the (110) directions below about 60 ev., with sharp focusing below 
about 15 ev. The presence of adjacent lines of atoms will enhance 

the focusing effect found in an isolated line of atoms. Although the 
hard sphere models used to demonstrate the phenomenon of energy 
focusing contain rather drastic assumptions, the same sort of effect is 
expected from models based on more realistic potentials. 


A channeling effect has also been saeomeneo 


This is the tendency 
for an energetic atom traveling in a (110) or (100) direction to be focused 
along a (110) or (100) line that is centered between surrounding (110) or 
(100) lines of lattice atoms. The energetic atom may travel a large dis- 
tance along a channel compared to the distance it would travel in some 
other direction. The rate of energy loss along the channel is highly 
dependent on the potential function which governs the interaction, the 
energy of the atom, and its deviation from the line of focusing. From 
geometrical considerations, one would expect (100) channels to dissipate 


more energy per unit length than would be dissipated by (110) channels. 
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Channels, energy chains, and the displacement chains, then, may 
produce vacancies which are separated by large distances from their 
corresponding interstitials and sputtered atoms. Thus, chains and 
channels are competing mechanisms which produce the same effects. 
Energy and displacement chains should not be significantly affected by 
the presence of vacancies. Vacancies in the lines of atoms surrounding 
a channel should tend to divert channeled energetic atoms into the 
corresponding chains. Interstitial atoms should tend to disrupt channels 
and chains. 

It appears that the chain and channel mechanisms are the dominant 
factors in the radiation damage event. The conditions for stability and 
the lengths of the chain/channels will then largely determine the final 


lattice configuration. 





3. Theoretical Work by Others 

Gibson, et al. a have studied radiation damage events by inte- 
gration of the classical equations of motion of the atoms in a small 
crystallite &1000 atoms). A central repulsive potential of the Born- 
Mayer type, A exp (-r/b), was assumed for the interactions between 
atoms. In order to contain the events it was necessary to restrict the 
energy to <400 ev. Calculations were made on a high speed digital 
computer. 


The computer program developed by Gibson, etal. 1 


keeps a 
record of position co-ordinates and velocity components for each atom 
in the crystallite. The orbits of the atoms are arbitrarily divided into 
segments with respect to time, i.e. time steps. In order to calculate 
the trajectories of the atoms for a particular time step, the atoms are 
considered sequentially. The resultant force on an atom under con- 
sideration is calculated from the position of the atom and the positions 
of the surrounding atoms. The atom is then moved for one time step 
under the influence of this resultant force. When the calculations are 
completed for all atoms for a given time step, calculations for the next 
time step are begun. 

If the repulsive force between adjacent atoms were the only force 
simulated, the crystallite would expand without limit. Therefore, it is 
necessary that an attractive force act on the surface atoms of the crystal- 
lite. If all collisions were elastic collisions, and if the original energy 
(for a 100 ev. run) were completely contained in the lattice, the mean 
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energy of the atoms would be about 0.1 ev. after thermalization. 
This corresponds to about 900° C. as the final temperature of a sys- 
tem which started at 0° C. Therefore , a viscous damping force was 
applied to surface atoms to dissipate energy. 

This model is useful in studying low and moderate energy (< 400 ev.) 
effects in radiation damage events. The effects studied were stability 
of lattice defects, anisotropicity of energy transfer (energy chains), 
and displacement chains. Interatomic potential and surface effects 
were also studied. 

In the Gibson model, energy chains play a dominant role in the 
dissipation of energy from the initial knock-on. Frenkel pairs are 
shown to be stable for distances of separation greater than 21/2-4 
times the nearest neighbor distance, depending upon the orientation 
of the pair with respect to the lattice. Interstitials were studied ex- 
tensively and found to be stable only in the split interstitial configura- 
tion. 

Oen, et ee have investigated the stopping of energetic atoms 
(1-100 kev.) by solids. The model is based on binary elastic collisions 
between a primary knock-on and successive target atoms which are 
initially at rest. Lattice effects are neglected; the target atoms are 
initially randomly distributed throughout the target material. 

The ultimate goal of the work by Oen, et al. is to discover the 
appropriate interatomic potential. The immediate aim is to calculate 
histories of many primaries which may be considered to simulate a 
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beam of incident atoms or ions. 

Trajectories and histories are calculated only for the primary 
knock-ons, i.e. any effect produced by sedondary knock-ons is 
neglected. The history of each primary consists of: 

1) R, the distance from the initial position to the final 
position of the primary, 
2) X, the component of R parallel to the initial direction of 
motion of the primary, 
3) P, the component of R in a plane perpendicular to the 
initial direction of motion of the primary, and 
4) L, the length of the path traversed by the primary. 
Each run consists of the calculation of histories for enough primaries, 
usually 1000, to insure that statistical fluctuations in the mean values 
of R, X, P, and L were acceptable (<€ 3%). 


The potential used was an exponentially screened Coulomb, or 


Bohr potential: Vcr) = Ep exp yA) 
6 a 
b 


where Eg, = CZ Z, C/A, and a,= kau Az +zZ*y* 

where ary = first Bohr radius (0.529 A). The parameter k may be adjusted 

to achieve agreement with experimental results. The hard sphere approxi- 

mation to the Bohr potential was also used in some cases for comparison. 
Results for the cases of 10.5 and 60 kev. Na24 ions incident upon 

aluminum targets were compared with the experimental results of Davies & 


Sims.” No serious discrepancy is found in the mean penetration, X. 
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However, the experimental distribution of penetration is more "skewed" 
toward lower penetration than the distribution of penetration calculated 


by Oen, etal. 
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4. The Problem 

The binary collision model has been the basis of all analytical 
studies of radiation damage and sputtering. Many machine calculations 
have been made in recent years using an N - body collision model. 
However, no serious attempt has been made to determine whether or 
not a two-body approximation would yield satisfactory results. 

Our problem, then, was to build a model based on the binary col- 
lision whose logic would be sophisticated enough to adequately repre- 
sent the radiation damage problem. The nature of the results obtained 
would then definitely establish the degree of confidence to be placed 


in any binary collision approximation models. 
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5. The Binary Collision 

With minor modifications, any repulsive potential may be used in 
the two body interaction. An eroded form of the potential is used, i.e. 
the potential has a definite radius of effect. We have taken this radius 
of effect (or sphere of influence) as an input parameter, SPHI. There is, 
then, no interaction between atoms whose Separation is greater than 
SPHI. 

All applications of the program to date have used an eroded Born- 
Mayer potential, V = A exp (B(l-r /SPHI)) for re SPHI. The constants 
A & B were taken from Gibson et al. for three different Born-Mayer 
potentials , termed simply Potentials, 1, 2, and 3. The constants for 


these three potentials are: 


POTENTIAL A (ev.) B 
] .0392 Logo 7 
2 .0510 13.00 
3 - 1004 10.34 


The interaction starts with the two atoms at SPHI distance of sepa- 
ration with all the energy as kinetic energy. All co-ordinates and velo- 
city components are transformed to the center of mass system, where 
the equations of motion are: 

1) dt = (2/m' (E-V-L2/2m'r2))7172 ar 


2 6c (L/m'r2) dt 


PS 


_- 





where t = time 

reduced mass 

kinetic energy 
potential energy 

total angular momentum 
distance of separation 
angle of deviation 
impact parameter 
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The equations are integrated from r = SPHI to r = CPA, the distance 
of closest approach. Since the interaction is symmetrical in the center 
of mass system, the definite integrals of the above equations are one- 
half the time of interaction and one-half the total angular deviation, 
respectively. 

CPA is found by an interation of the formula: °° 
(cpa)? (1 - V(CPAYE,) - g7=0 
The first approximation to CPA is taken as SPHI. A four point Gauss- 


| 
Legendre quadrature is employed to evaluate the integrals. “? 
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6. Development of the Model 

When the possibility of a feasible two body approximation to the 
radiation damage problem is considered, two major problems are im- 
mediately encountered. 

The atoms would, of necessity, be considered in some sort of 
ordered sequence, but all calculations for a single atom could not be 
carried through before the other atoms are considered since this would 
obviously leave out of the calculations the effect of collisions suffered 


by other atoms. This first major problem led to the time step approach, 


and eventually to the use of an absolute time of an atom's last collision. 


The second major problem is to develop a mechanism which will 
handle simultaneous collisions within a two-body approximation. 

After the program was developed, we discovered that the lattice 
would not hold together which necessitated a further modification. 
This difficulty is only partially resolved. 

The first question considered in the development of the model was 
to mathematically create a face centered cubic lattice in table form, 
together with a method which tags each atom so that its subsequent 
movements and initial position could be found. Figure 2, APPENDIX II, 
shows the ordered arrangement found in a (100) plane of a fcc lattice. 
The distance between two atoms along the (100) direction is taken as 
two units. It will be noted that if the total X dimension, IX, is odd, 
then there are an equal number of atoms in the line y = 0 and the line 
y= 1. By specifying IX odd, the number of atoms in a line y = constant 
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is known to be (IX +1)/2. This value, NLINE, can then be multiplied 
by the total Y dimension plus one (IY+1), and the number of atoms in 
the whole plane, NPLANE, is obtained. It will be noted from Figure 3 
that the next plane, one Z unit removed, is akin to a photographic 
negative of the plane under consideration. However, the ideas of 
NLINE and NPLANE are still valid. The total number of atoms in the 
lattice can then be found by multiplying NPLANE by the Z dimension 
plus one (IZ+1). If the atoms are numbered sequentially as shown in 
APPENDIX II, the initial co-ordinates of their sites can also be found 
by the method shown in APPENDIX II and APPENDICES IV & V (BOX 31). 

The process outlined above will create and number a 1500 atom 
lattice in approximately one second. An advantage of this method is 
that a rectangular lattice of any size or shape desired can be created 
by specifying IX, IY, and IZ. The only restriction is that IX must be 
odd. 

The atoms are now considered sequentially, for a period of one 
Time Step Length (TSL), and all collisions that will take place before 
the end of the time step are allowed to happen. All atoms are then 
moved to points in space and time corresponding to the end of the time 
step, and the process is repeated for the next time step. 

The length of a time step was originally calculated so that the 
most energetic atom in the lattice could move one lattice unit (1.807 
Angstroms for copper) in one time step. This has been modified so 
that it is now able to move any desired fraction, TFAC, of one lattice 
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unit. TFAC is, naturally, one of the input parameters of the problem. 
Since the energy of the most energetic particle, ENERGY, in the 
lattice decreases during the course of the run, the TSL should be re- 
calculated periodically for greater efficiency. Provisions were made 
to recalculate the Time Step Length every JPB (an input parameter) time 
steps. The TSL is also recalculated every time another particle is shot 
into the lattice, since it is assumed that this particle will be the most 
energetic one present, and it was not thought feasible to wait for the 
next regular recalculation. A table of the first 30 Time Step Lengths, 


and the time step of recalculation, is kept for reference. 


As now constituted, there are one external and two internal methods 


used to determine the number of time steps which complete arun. The 
external method is the only one available to the operator once a com- 
putational run has commenced. 

First, a variable limit is set on the maximum number of time steps 
(MNTS, an input parameter) that may be taken. The calculation must 
cease when this number of time steps has been completed, and the 
results are then obtained. 

Second, the program, at the beginning of each time step, decides 
whether or not it should continue for another time step. The decision 
to continue is made if the energy of the most energetic particle (called 
ENERGY) in the lattice (excluding all those that have left the lattice) 
is greater than a certain threshold energy, THERMAL (an input parameter 
explained below). If ENERGY is leg$ than THERMAL, the run is 
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considered completed and the informational output is given. 

Third, the operators may decide to shut down operation before 
either of the two internal limits have been reached. This may be done 
at any time, simply by throwing a switch, and the computation will 
cease at the beginning of the next time step. 

The external mode of problem completion should not be the normal 
operating condition, nor should the number of time steps be allowed to 
approach MNTS under most conditions, since either method still leaves 
atoms in the lattice with more energy than the quantity THERMAL. 
Usually the calculation should cease of its own accord, with the 
second method given above as the basis for the internal decision. 
Other factors, such as running time, mechanical failures, etc., may 
necessitate the use of the external method, but these are to be avoided 
if possible. 

A method of program regeneration was also devised so that the run 
could be discontinued at the start of any time step and then restarted at 
that point at some later time. This has definite value where continuous 
periods of available machine time are not long, or when mechanical or 
electrical failures destroy the calculations. The information needed for 
regeneration can be obtained at the beginning of any time step and the 
calculation continued (or stopped if desired). 

In this manner, if the calculation has been running for three hours 
before a failure of some kind occurs, the only time that need be re- 
peated on the machine is the time elapsed since the last regeneration 
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information was put on magnetic tape. It should be noted that the re- 
generation need not be undertaken at once, but may be delayed until 
convenient. This regeneration process has great practical importance 
to the user of the program, but is of no theoretical value. 

Many atoms are given only a cursory examination during a time 
step, while the more energetic atoms are subjected to a highly detailed 
treatment. The use of two energy thresholds in the model simplifies 
this choice for a particular atom. 

ETHRESH, the upper energy threshold, is the energy required for 
an atom to move from its site to a site near it, i.e. the displacement 
energy. Calculations by others!” have shown this displacement 
threshold to exist, and also that it is in the neighborhood of 20 


electron volts. They have also shown! 


that this displacement thresh- 
old depends upon direction, and is lower in the (100) and (110) direc- 
tions than in other directions in the lattice. This directional depend- 
ence of the displacement threshold has not been built into the present 
model, but it should be considered in any future development of the 
model. 

THERMAL, the lower energy threshold, is an arbitrary dividing 
line imposed upon the model, below which atoms are quite arbitrarily 
placed on sites or formed into interstitial pairs with other atoms and 


not allowed to move. Atoms above THERMAL are allowed to participate 


normally in collisions. 
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There are, then, three distinct energy classes among the atoms in 
the lattice: 


1) Those with energy above ETHRESH, 





2) Those with energy between ETHRESH and THERMAL, and 
3) Those with energy below THERMAL. 
The three classes will be discussed separately. 

The atoms above ETHRESH do not occupy a site in the lattice, 
collide with other atoms as necessary, and obtain information about 
their subsequent motion only from the interaction subroutine. 

The atoms between ETHRESH and THERMAL are forced to occupy 
specific lattice sites, but are allowed to vibrate around these sites. 
They cannot move to a new site because their energy is less than the 
displacement energy. These atoms are allowed to move undisturbed 
until they have a collision, and then after the collision has been com- 
pleted, if they lose a specific percentage, TURN (an input parameter) , 
of their original energy their velocity vectors are changed to point back 
toward the site occupied by the atom. 

We found prior to the use of this turn-around method that the lattice 
had no inherent cohesiveness. This method holds the lattice together. 

To speed up computation time, atoms with energy below THERMAL 
are forced to occupy a fixed site, or to share an interstitial site with a 
partner, and are not allowed to move. This is a reasonable and valuable 


approximation for atoms of low energy, since they will not have an 
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appreciable effect om the re of the computation. The energy and 
associated velocity components of these thermalized atoms are kept 
intact for use in possible future collisions. 

After the decision is made to continue for another time step, and the 
Time Step Length has been recalculated, the atoms are considered sequen- 
tially, The information associated with the atom under consideration 
(primary atom) which is stored in the LB array is first inspected. If the 


atom has: 


Res 
{ 


1) An energy mh to or less than THERMAL, or has 
2) Been through the current time step, or has 
3) -Left the lattice, 

we move on to consideration of the next atom. 

We will account for all collisions if we consider only those atoms 
with energy greater dias THERMAL, since atoms with energy less than 
THERMAL are not allowed to move unless hit. Atoms that have already 
been through the time step have completed all interactions which occur 
during the time step and should not become involved again until the 
next time step. Atoms that have left the lattice will not be considered 
for the remainder of the problem. 

A primary atom with energy above THERMAL is finally selected, or 
the computation shuts down. A list of all other atoms it is possible for 
this atom to hit in this time step is required. To shorten the computa- 
tion, a mathematical box is constructed around the atom in question 
(hereafter referred to as atom M) and a list (the NUM list) is made of 
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all atoms whose positions (centers) are inside the box. If the box is 
made large enough, it will be impossible for atom M to hit any atoms 
outside the box, so we need consider only those on the NUM list. At 
present, the box is given sides of 3.02 units, or + 1.51 units. Initially, 
the sides were + 1.01, but this has been found to be too small, and 
obvious collisions were missed. Atoms that have been through the time 
step are not included on this list, since their co-ordinates are relative 
to the end of the time step rather than the start of the time step. Atom M 
is placed first on the list for programming facility and easy reference. 

The impact parameter between every other atom on the NUM list and 
atom M is then calculated, these distances comprise the DSTANCE list. 
The time at which each atom reaches this distance of closest approach, 
and the time at which the atom and atom M are SPHI (sphere of influence 
of an atom) distance apart are also calculated. These times are placed 
on the TMIN and T lists, respectively. 

During the course of the calculation of the DSTANCE, TMIN, and T 
lists, several tests are performed on these quantities. If the DSTANCE 
for an atom is larger than SPHI, no collision will occur. If the relative 
velocity is less than 1076, then the atoms are moving very slowly 
relative to one another, and the collision is ignored. Several tests are 
then performed on each T. If the time T is less than the start of the pre- 
ceeding time step, it would have been considered in this previous time 
step, and is therefore treated as a spurious collision. If the time T is 
greater than the end of the current time step, it will be re-examined in 
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the next time step. The times of the last collision of atom M and the 
other atom are then checked, and if the time T, on an absolute scale, 

is less than either of these, the collision is spurious since it would be 
starting before one that had already taken place. The last test performed 
is whether the time T is within .00001 (TSL) of the end of the time step. 
If so, we feel that it would be better to consider this collision in the 
next time step, since to consider it now would mean that the end of the 
collision was far into the next time step, which is unsatisfactory from 

a practical viewpoint. 

If, after this series of eliminations, there are no atoms left, then 
atom M is moved to the end of the time step, and the next atom in se- 
quence becomes atom M. 

If atoms remain on the list we select the ones that atom M will hit. 
Two alternative methods of selection were considered by the authors: 

1) Select the one, or ones, that are closest {(i.e. a space- 
wise selection), or 
2) Select those that will collide with atom M first (i.e. a 
time-wise selection). 
Only the second method has been used. 

The space-wise selection would surely include all the hardest 
collisions, and would probably be easier from a programming viewpoint. 
However, a time-wise selection insures that atom M will hit first what 
it really should hit first. Atom M may thereby miss a hard collision 
with another atom because of another previous soft one, but we 
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currently feel that this is a truer picture of what actually happens in the 
lattice during the radiation damage process. We feel that both methods 
should be tried in future development of the model. A study as to the 
relative feasibility of both methods would enhance the usefulness of the 
model. 

Atoms are selected that have the smallest times, T. Since we assume 
that two atoms have no interaction if they are separated by a distance 
greater than SPHI, the time T is the actual time of collision. We feel 
that some leeway should be available here, therefore all atoms whose 
times T are within a certain amount (CUTOFF) of a time step from the 
minimum times are also selected. These atoms are placed on the KHIT 
list, with atom M first on the list. Since CUTOFF is an input parameter, 
adjustment to a satisfactory value may be accomplished easily. Use of 
the CUTOFF parameter also eliminates some of the objections to the 
time-wise selection process. 

The initial velocity components and energies of all atoms on the 
KHIT list are then stored in the SAVE array because they will be needed 
for the velocity scaling processes. 

It must be emphasized that atom M hits all those on the KHIT list 
almost simultaneously. This necessitates some sort of energy and 
velocity scaling process, or a general solution to the N body problem. 

We feel that if the two body approximation is to have any validity 
atom M must enter into each of the two body collisions with its initial 
energy, but the program must scale the results in some fashion. Assume, 
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for example, that M initially has an energy of 100 ev. and hits two 
atoms simultaneously with identical impact parameters. It is more 
feasible to assume that M hits each one with an energy of 100 ev. 

than to assume that M hits one with 100 ev. and hits the other with 
whatever energy is left, say 37 ev. 

Naturally, with M hitting each atom separately with 100 ev. the 
struck atoms will receive too much energy, we must therefore scale the 
results. 

In the actual collision process of the model, then, atom M and the 
next atom on the KHIT list are moved to their calculated positions at the 
time T for the atom being hit. The collision is then calculated and the 
struck atom is given its new position co-ordinates, velocity components, 
and energy. Atom M and the next atom on the KHIT list are then moved 
to the positions corresponding to their actual time of collision, T, and so 
forth. This is repeated until atom M has "struck" all the atoms on the 
KHIT list. 

To conserve both energy and momentum in any scaling process 
would constitute a general solution to the N body problem. We, there- 
fore, elected to conserve energy rather than momentum, feeling that 
this would give us a more realistic insight into the radiation damage 
problem. 

There are four general scaling cases that must now be considered: 

1) Atom M "lost" more energy than it started with, 
2) Atom M "lost" exactly all the energy it started with, 
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3) Atom M started with an energy apove ETHRESH and did not 
"lose" all of its energy, and 
4) Atom M started with an energy below ETHRESH, but above 
THERMAL, and did not lose ail of its energy. 

These cases will be discussed in order. 
CASE I 

Atom M is allowed to lose al! its initial energy, and its velocity 
components and energy are set equal to zero. The energies and 
velocities of the struck atoms are then scaled to conserve energy. The 
energy scaling factor used in this case is simply the total energy (of 
those on the KHIT list) before the collisions divided by the total energy 
after the collisions. The velocity scaling factor is the square root of 
the energy scaling factor. The energies and velocities of all atoms on 
the KHIT list, except atom M, are scaled hy these two scaling factors. 

Both energy and the absolute velocity directions are conserved by 
this process. 
GaSEnI 

Atom M is allowed to lose all of its initlal energy, and its velocity 
components and energy are set equal to zero. Here, however, energy 
scaling is not necessary since energy has already >een conserved. 
CASE Ill 

Here it is not necessary to scale the energies of the struck atoms, 
since atom M did not lose all of its initial energy, so the struck atoms 
retain the energies and velocities they received in the two body collision 
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process. The energy of M is found by subtracting the energy gained by 
the struck atoms from M's initial energy. M's velocity components 
must now be found. 

The sums of the individual velocity components calculated for M, 
but not given to M, are taken. Call these sums S);, S2, and 53 for the X, 
Y, and Z directions. Let R be is? + sé + aayl/?, Tne magnitude of M's 
final velocity, V, is known since the energy is known. The X velocity 
is then taken as (V/R) (S ) , and similarly the Y and Z velocities are 
(Vv /R) (S.) and (V/R)(S3) respectively . 

CASE IV 

Two sub-cases must be considered. First, if atom M does not lose 
TURN % of its initial energy, the treatment used in CASE III will be 
followed. 

Second, if atom M does lose a significant portion of its energy, its 
velocity components are changed so that the atom is directed back toward 
the site it occupies. The energy of M is again its initial energy minus 
the energy gained by the struck atoms. The scale factor used is 
(Velocity/distance from site). This is multiplied by AX, AY, and 
JA Z to give the velocity components. AX is simply the X co-ordinate 
of the site minus the X co-ordinate of M. 

In all cases, atom M's spatial co-ordinates are found by averaging 
the individual co-ordinates calculated for M in the X, Y, and Z directions 


(i.e. dividing the sum by NHIT). 
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If atom M is a primary atom (i.e. one of those considered sequen- 
tially for this time step), then the KHIT list is duplicated by the LAST 
list. Atom M will then proceed to go through its designated collisions, 
look for more collisions, and go through these until it completes all 
collisions in this time step. 

When atom M has completed all its possible collisions for this time 
step, each atom that M hit (found in order on the LAST list) undergoes 
the same process and attempts to hit other atoms. The atoms hit by M, 
or others on the LAST list, are not inspected for possible further colli- 
sions in this time step. We feel that these "secondary" atoms would 
not have enough time to have a collision in this time step, and if they 
did, the collisions will be calculated in the next time step. At this point, 
a possibility does exist that some reasonable collisions are excluded 
by the process, but extensive calculations have not revealed any 
examples of this exclusion to date. 

The absolute time of an atom's last collision is also calculated 
during the collision processes described above. For each of the struck 
atoms this time is the total time from the start of the computation to the 
start of this time step plus the time, T, of collision (T is measured from 
the start of the current time step). 

The time of M's last collision is taken as the start of its collision 
with the last atom on the KHIT list (computed as above) plus a certain 
percentage (PERCENT, an input parameter) of the present time step 
length. M is thus prevented from participating in further collisions 
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for a definite period of time after the start of its last collision. We 
feel that this approach is a valid one, however study of this aspect of 
the model should be continued to determine its true usefulness. 

Sites must now be found for those atoms on the KHIT list whose 
pre-collision energies were above ETHRESH and whose energies are 
now below ETHRESH. 

If the rounded-off co-ordinates of an atom in this class (call it J) 
give the co-ordinates of an actual lattice site, then the atom will 
occupy that site, possibly as part of an interstitial pair if the site is 
already occupied. Atom J will then be allowed to vibrate around its 
site until its energy falls below the THERMAL cutoff. If J's rounded-off 
co-ordinates are not the co-ordinates of a site, then the nearest neighbor 
sites are found and the site closest, space-wise, to atom J is chosen as 
J's vibrational center. Again, this may result in the formation of an 
interstitial pair. Since one of the atoms has an appreciable energy, 
however, this interstitial pair should be destroyed in the next time step. 

The site decided upon may already be occupied by an interstitial 
pair: if so, a "triple interstitial" is formed. An error counter is then 
increased (MISTAKE (5)). The possibility of "triple interstitials" could 
be eliminated by more extensive programming, but an evaluation should 
first be made of the seriousness of this error. The frequency of forma- 
tion of these "triple interstitials" has been small thus far, not more 
than two or three per run, which is tolerable. If this frequency reaches 
an intolerable level in future development, the model must be modified 
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to eliminate any possibility of their occurrence. 

All atoms on the KHIT list whose post-collision and scaling energies 
are below THERMAL are placed on the LATER list at this time. Use of 
the LATER list is discussed below. 

If one of the struck atoms on the KHIT list was a member of an 
interstitial pair, the other atom must now be considered. If this other 
atom was not hit, and the struck atom now has an energy above ETHRESH 
(a very unlikely occurrence), then the other atom is placed back on the 
site, rather than being left in the split interstitial position. We feel 
that moving the atom back onto the site is more realistic than the 
definite error of leaving it in the split interstitial position. 

All calculations in the model are made on the assumption that the 
atoms are both space-wise and time-wise at the start of the current time 
step, except those that have completed the current time step. Therefore, 
all the atoms on the KHIT list are now moved back in time and space on 
their post-collision tracks to the start of this time step. 

The entire collision process, including the construction of a new 
NUM list is now repeated for atom M. This repetition will continue for 
M until M is unable, time-wise, to have further collisions in this time 
step. The next atom on the LAST list will then become atom M and the 
collision process will be repeated until this M is unable to have further 
collisions in this time step. This overall procedure is then repeated 


until the LAST list has been exhausted. 
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All the atoms on the LAST list, except those on the LATER list, 

(the original M is the first atom on the LAST list) are now moved, space 
and time-wise, to the end of the time step. These atoms have now com- 
pleted the time step and will not enter into the subsequent calculations 
until the next time step. 

The LATER list is now utilized to place the atoms that have partici- 
pated in this sequence of collisions and whose energies are below 
THERMAL on sites, either as members of interstitial pairs or by them- 
selves on previously vacant sites. 

The co-ordinates of the atoms on the LAST list are now examined. 

If an atom has left the lattice, the time step number and the particular 
face (see APPENDIX I, NFACE) of exit is recorded in the LB array. 

This completes the consideration of atom M (the original M) for this 
time step. The next atom (M + 1) is then considered in the same manner 
until all atoms have been through the procedure (and, consequently, the 
time step). 

After all the atoms have been through these processes, all un- 
affected atoms are moved to the end of the time step, which completes 
the time step. A new time step is then started, or the run is shut down 
by one of the three processes mentioned earlier. 

A numerical example of the order in which atoms are considered may 
be of some value here. Assume atom #73 (M) hits four atoms simulta- 
neously, atoms 75, 79, 87, and 65. Assume further that #73 hits them 
in the order: #79, #65, #87, #75. Atom #73 will hit these four, then 
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repeat the process for further collisions. Assume #73 hits #80 and #67 
next, but has no further collisions after that. Then, atom #79 will 
become the next M since it is second on the LAST list, #65 the next M, 
and after #75 has completed its possible collisions, if any, the regular 
sequence will be taken up once more. Atom #73 was the last one con- 
sidered sequentially, therefore atom #74 now becomes the primary atom, 
atom M. 

As mentioned above, atoms whose energies fall below THERMAL are 
either placed alone on vacant sites, or formed into split interstitial pairs 
with other atoms near them. This is done only after a primary atin , to- 
gether with thoseonits LAST list, has been completely considered. If 
this were done at the end of every collision process, it would severely 
hamper the multiple collision process, and impose a great deal of strain 
on the credibility of the argument behind it. 

When a list of the nearest neighbor sites is formed, if the atom in 
question (call it MJ) is within one unit of a face of the crystallite, some 
of the geometrically nearest neighbors normally found will not be present. 
There are 22 special cases for each one of the two general cases (MJ's 
co-ordinates are a site or non-site). Sixteen of these special cases 
for each concern the edges and corners of the crystal. In our model, 
the 32 (total) possible cases for atom MJ near an edge or corner are not 
considered separately, but are included in the 12 special cases of when 
MJ is near a lattice face. When one of the 32 unconsidered special 
cases does arise, an error counter is increased. Inspection of all 


32 








Et ne ee 
wwe i 
CF Remit my amie lh bag 
lt ae 1) 
el cece lll hdl 
| eh a a mv cE 
a af ewe’ of etl: Qi s Gp ag 
ed ee & Gee 6 ee a ep 
yee TN ae AP Eee 
5) ee nee | ee ih dimmed bocca Fs eee» 
—_—_———— eae” gerne eee Oe teemmiveeme 









— « =? & —— i ee a 
= == Gas « —— © 1 mh eee gem ane Apel ee © 
—— - Rw eee ee 





Ee — tail — tee “Sent? Be 
= — eae =6—8 ee 


» 


preliminary results to date have, however, shown NO discrepancies as 
a result of this simplification at the edges and corners of the lattice. 

After the list of nearest neighbor sites, the MOON list, is formed, 
these sites are inspected for possible vacancies. If only one vacancy 
exists, atom MJ is placed on that site and the appropriate corrections 
are made to the LB array. If more than one vacancy exists among MJ's 
nearest neighbor positions, then one of these must be selected for MJ's 
occupancy. We feel that this selection should be made on a space-wise 
condition of proximity, rather than on any time-wise condition. There- 
fore, atom MJ is allowed to occupy the closest vacant site. If there is 
more than one vacancy at this smallest distance, MJ will be placed on 
its own site, if available, or the selection will be made in a random 
manner. The choice of a particular vacancy in the case of multiple 
vacancies is somewhat unsophisticated, however at this stage of devel- 
opment in the radiation damage field, we feel that no process of a more 
sophisticated nature can be justified. 

If no vacancies exist among atom MJ's nearest neighbor sites, then 
a prospective interstitial site and partner are chosen. If atom MJ has 
no energy or velocity, this decision is made on a space-wise basis, 
i.e. the closest site is chosen. If atom MJ does have energy, then the 
selection is made on a time-wise basis, i.e. the site that atom MJ 
passes first. 

Once the selection of a prospective interstitial partner has been 
completed, the atom (call it MM) which occupies the site is found. 
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Atom MM's nearest neighbor sites are then found, and the list is 
examined for vacancies. If any vacancies exist, then atom MM is 
moved to one of these vacant sites by the same method discussed 
above for MJ. Atom MJ is then moved onto atom MM's recently vacated 
site and no interstitial is formed. 

If, however, atom MM finds no vacancies, then an interstitial pair 
must be formed. We feel that of the three types of interstitial pairs 
mentioned earlier, the split interstitial is the most plausible. Indeed, 
Gibson, et al. have shown the split interstitial to be the only stable 
configuration. !? The AX, AY, and AZ distances between atom MJ 
and the site atom MM occupies are found. The largest of these is 
taken as the axis of the interstitial pair. Each atom is then placed on 
its split position, 0.5 — the site along the axis. Once the 
appropriate corrections are made to the LB array, the formation of the 


pair is completed. 
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7. Operating Procedures 

The model consists of three separate but inter-connected programs, 
MASTER, SLAVE, and RON. Machine memory limitations forced the devel- 
opment of three, rather than one, programs. 

The main program, MASTER, is the computational program and does 
all of the calculations necessary to complete the run. The two second- 
ary programs, SLAVE and RON, merely use the binary output of program 
MASTER as input and put the results in a semi-analyzed readily legible 
form. The orbits (or intermediate positions of the moving atoms) are 
given by program RON while program SLAVE tabulates and categorizes 
the final positions, energies, etc. of the lattice atoms. Again, be- 
cause of memory limitations, the two secondary programs may not be 
combined. 

If a binary version of program MASTER is obtained with the FORTBIN 
compiler, the storage required may be found by listing of the binary tape. 
The maximum lattice may then be computed. The FORTMAP compiler 
should be used to compile the large version of MASTER since it requires 
less memory than the other compilers. 

Once program MASTER is compiled, a run may be started. MASTER is 
called into memory and the run begins. MASTER calls for BCD input on 
Tape 2, and this is normally transferred to the card reader. BCD output 
is given on Tapes 3 & 4. Tape 4 is normally transferred to the typewriter 
for output (see BOX 10, APPENDIX V). The BCD output on Tape 3 is 
temporary and may be eliminated when the program is working correctly. 
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This will also make more memory space available. 

Program MASTER writes binary output on Tapes 6, 7, and 8. The 
MASTER output on Tape 6 is the binary input to program RON. The MASTER 
output on Tape 7 is the regeneration output and is used only by the MASTER 
program. The MASTER output on Tape 8 is the binary input to program 
SLAVE. MASTER will automatically rewind each of the three binary tapes 
when output for a particular tape is completed. 

When MASTER has completed computation program SLAVE should be 
called from our system library. SLAVE calls for BCD input (comment 
cards, etc.) on Tape 2, and this is normally transferred to the card 
reader. SLAVE will then read Tape 8, compute, and write its BCD out- 
put on the designated tape. Program SLAVE, at this time, has PRINT 
statements for all output except that mentioned in BOX 1 of the SLAVE 
part of APPENDIX V, and therefore the output will be on Tape 4 unless 
otherwise directed by the operator. 

When SLAVE has completed its run, program RON should be called 
from the system library. The input and output designations for RON are 
the same as those for SLAVE except that RON reads Tape 6 rather than 
Tape 8. 

If, for some reason, MASTER does not complete its run, and Tape 6 
is not rewound, program RON must be run with Jump Switch 1 in the set 
position. This will insure that the correct END designation (in this case 
a -1) is placed at the end of the information on Tape 6 and that the tape 
is rewound. 


36 










—_ <3 4 «= oe 
—— 
2S Ye 
——_ << «= a alee | 
ae 668 eee Gees ep 
—  —— «a <a? oe 








—- «© —_ 2 Ee, 


<<: f° + |, | 





(— ee — ee a2 eS 15, 
> - acta a eta, aisle iiss se 
a 7 Fe 
; «* — — at fa Gp Gee € fF 


aE 





8. Parameter Adjustment 

The two input parameters, ETHRESH and TURN primarily govern the 
cohesiveness or solidity of the lattice. At ETHRESH = 3 ev., for 
instance, some 400-500 replacements occur with a 100 ev, shot. With 
ETHRESH raised to 15-25 ev. the number of replacements drops to con- 
siderably less than ten, depending on the location and direction of the 
initial knock-on atom. The best working values seem to lie in the 
15-25 ev. range, as preliminary results in this range show good correla- 
tion with the work of others. More accurate values of ETHRESH and 
TURN should be obtained in the future from concurrent N-body calcula- 
tions for use in the model. 

Adjustment of the TURN parameter also governs the length of energy 
chains. If TURN = 0.0, the energy chains are almost non-existent. 
Atoms vibrating on their sites are then re-directed toward their sites 
after their first collision. Since this first collision is usually a slight, 
glancing one, very little energy is transferred and the energy chain is 
not seen. If, however, the TURN parameter is set too high, displace- 
ments occur among atoms with energy< ETHRESH. Investigations thus 
far seem to indicate that a value of about 0.1 is the optimum value for 
TURN. 

THERMAL has little effect on the larger aspects of radiation damage 
i.e. chains, channels, displacements, ring displacement, etc., but 
has a large effect on machine running time. By changing the value of 
THERMAL from 0.1 ev. to 0.5 ev., the machine time can be reduced by 
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75% or more. The very low (below THERMAL) energy vibrations are no 
longer seen, to be sure, but since this is not one of the areas of pri- 
mary interest, the advantage of drastically decreased machine time far 
outweighs the inherent disadvantages. 

When the penetration of higher energy (above a few hundred ev.) 
initial knock-ons is the primary effect under consideration, THERMAL 
may be set as high as 3-5 ev. For initial knock-on energies of 100 ev. 
and below, THERMAL should be 0.1 - 0.5 ev. to achieve a reasonable 
degree of thermalization. 

Investigations have shown that the radius of effect (sphere of in- 
fluence), SPHI, should lie within the range 2.5-2.6 Angstroms (for 
copper). Values much outside this range (2.2 say) give rather meaning- 
less results. We feel, as a result of the studies we have made, that 
the best value is probably slightly less than the nearest neighbor dis- 
tance. For copper, the nearest neighbor distance is (2) 1/2 (1.807). 
We have used the value (1.414) (1.807) when not specifically investi- 
gating SPHI effects. We feel that any value larger than the nearest 
neighbor distance will produce invalid results. 

Studies of the optimum values for PERCENT and CUTOFF are not 
completed at this time. Preliminary investigations show that values 
above 0.1 for PERCENT and below 0.1 for CUTOFF prevent significant 
collisions from taking place. 

The value of TFAC should be kept below 1.0, unless the normal 
value of 1.0 is used. A more accurate tracing of individual particle 
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movement is possible with smaller TFAC's, but machine time increases 
proportionately. A compromise must then be made between the running 
time available and the individual tracing accuracy desired. 

FMASS may be adjusted for the particular metal being used. At 
present the model is restricted to metals with face centered cubic struc- 
ture. The Born-Mayer potentials used are restricted to copper only 
since values of the interaction input parameters (the constants A & B 


mentioned above) are not available for other materials. 
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9. Preliminary Results 


Preliminary results to date show that the existence of channels, 
energy and displacement chains, vacancies, interstitials, and replace- 


ments depend quite heavily, and quite critically, on the values of the 


input parameters used. Before very satisfactory results can be obtained, 


more reliable values of these parameters must be determined. 

A run corresponding to Gibson, etal's. 17 Run No. 12 has been made 
with a slightly modified version of the program, using a 40 ev. primary 
knock~on. The primary knock-on was shot directly at one of the atoms 
on the front face (Gibson et al. start their primary knock-on inside the 
crystallite), at an angle of 15° away from the normal (100) direction. 
The values of the various parameters for the run were: 

DIX xXIY xIZ=11R 1526 

2) ETHRESH = 20.0 ev. 

3) THERMAL = 0.5 ev. 

4) CUTOFF = PERCENT = .15 (i.e. 15%) 
5) SPHI = 2.5 


6) TFAC = 1.0 


7) TURN = 0.1 
Our results (see Figure 1) are similar to those obtained by Gibson, 
etal. in that the energy chains play a dominant role in the dissipation 
of energy. The (100) chain of our run is not as prominent as Gibson, 
et al's., but strong focusing is observed in the (110) directions. Our 
(100) chain was terminated prematurely by the velocity turn-around 
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Figure 1. 
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mechanism employed for lattice cohesiveness. 
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10. Recommendations 

Many ideas and plans occurred to us during the development of the 
model to its present level that we have not been able to institute be- 
cause of time limitations. These unrealized plans and ideas are offered 
as possible guidelines for future development. 

We are convinced at this time that the programming errors in the 
model are concerned almost wholly with the LB information array. As we 
did not anticipate the great dependence upon the LB array when we began 
the problem, we feel that the LB aspects of the entire model (all three 
programs) must be overhauled, and more orderly storage of the informa- 
tion should be investigated. 

Once this overhaul and possible re-arrangement of the LB array has 
been accomplished and the model is working as designed, a more com- 
plete study should be made of the effects of parameter adjustment. 

At present when an atom on the LATER list doesn't find a vacancy, 

a prospective interstitial partner is chosen from the MOON list. A 
search is then made for vacancies near this partner. If none are found, 
the interstitial pair is formed. We strongly recommend that when the 
prospective partner finds no vacancies, the rest of the atoms occupying 
nearest neighbor sites should also be inspected for vacancies near them. 
If, again, no vacancies are found then the interstitial pair should be 
formed between atom MJ and its original prospective partner. 

To facilitate this, we recommend that the INVAC subroutine be 
split into two or more subroutines. There should be a separate 
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subroutine whose only function is to find the nearest neighbor sites, 
not necessarily by the present system. 

In the present model, if a vacancy occurs near an interstitial pair 
after the interstitial pair has been formed, this Frenkel pair will not re- 
combine unless the interstitial pair is hit by another atom. Provisions 
should therefore be made for periodic inspections of all interstitial 
pairs and any possible vacancies near them. A procedure very similar 
to that discussed above for making interstitial pairs could be used. 

We also feel that subroutines and functions should be incorporated 
to a much greater extent than at present, and that wider use should be 
made of comment cards in the programs. 

A study of the relative merit of a space-wise versus a time-wise 
selection process for the KHIT list has been suggested above and should 
be made. 

The ability to use atoms of two different masses should be built 
into the model. The actual use of two different masses with a Born- 
Mayer potential must be deiayed until appropriate constants for the 
potential are found. 

ETHRESH should be made a directional dependent variable rather 
than a constant as at present. Since a function able to represent this 
directional dependence of the displacement energy is not available, a 
table of thresholds associated with direction seems the most promising 


course of action. 
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Two memory saving ideas have occurred to us: 
1) The possible elimination of the TLAST concept, and 
2) The possible use of a method, already developed, which 
packs the co-ordinates and velocity components of an 
atom into three memory cells rather than the six required 
at present. This sub-program, SQUEEZE, developed by 
one of the authors (RLK), saves space at the expense of 
running time. 
The use of these two ideas would save considerable memory space, and 
enable the users to approximately double the size of the available 
lattice. 

The energy and velocity scaling processes used after the interaction 
(see BOX 20, APPENDICES III, IV, & V) are a first approximation to the 
problem and should be considered in that light. We feel that at present 
this is one of the weakest areas of the model. Considerable effort should 


be devoted to this section in the future. 
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11. Conclusions 

The model is not fully developed at this time, but we feel that 
some valid conclusions may be drawn from the results obtained thus far. 
The MASTER program (see APPENDIX III), as it exists now, still contains 
some minor programming errors but we feel that the logic is basically 
sound. 

The outstanding attribute of the model at this time seems to be the 
greatly decreased machine time necessary to achieve results compared 
to the machine time necessary for an N -body model. This decreased 
machine time makes it possible to use a lattice as large as the machine 
memory can accommodate. 

We do not feel that this model is sufficiently sophisticated at this 
time to state definitely the degree of confidence that should be placed 
in it. We do feel, however, that the preliminary results show that with 
more work the model can be made into a reasonable approximation to the 
radiation damage problem. We anticipate that our specific recommenda- 
tions, if carried out, will lead to a program which approaches this 


required degree of sophistication. 
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APPENDIX I 


SPECIFICATION OF UNITS AND ANNOTATED LIST OF VARIABLES 


The units used in the model for energy, mass, length, time, and 


velocity are given below: 


ENERGY 


MASS 


LENGTH 


TIME 


VELOCITY 


The electron volt (ev.), equal to 1.60207 x 107! 9joules é 


The Atomic Mass Unit (AMU), equal to 1.65983 x 10727 
kilograms. 


The Angstrom unit, equal to 10° *° tate! 
The jiffy, defined as 10714 seconds. 


Derived. One Angstrom/Jiffy = 104 meters/second. 


The variables given below are essentially in the order of their 


appearance in the program listing. However, certain basic variables 


are given at the beginning of the list. An asterisk (*) denotes an input 


variable. 
VARIABLE 


Beir, v) 


USAGE 
An array containing the position co-ordinates, velocity 
components, and energies of the atoms in the lattice. 
Lis the number of the atom. For I there are seven cases: 
1 The X co-ordinate in Angstroms 
2 The Y co-ordinate in Angstroms 
3 The Z co-ordinate in Angstroms 
4 The X velocity component in Angstroms/jiffy 
5 The Y velocity component in Angstroms/jiffy 
6 The Z velocity component in Angstroms/jiffy 
7 The energy in electron volts 
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VARIABLE 


LB (L) 


USAGE 


An array used for information storage, supplementary to 


the B array. Lis the number of the atom. Stored from 


right to left we find the following information: 


OCTAL 

POSITION 

i Bit 1 
Bit 2 
Bit 3 

2-5 

6-9 

On de 

14 


MEANING 
1 if the atom has completed the current time 
step. 
0 if the atom has not completed the current 
time step. 
1 if the atom has an energy > THERMAL. 
0 if the atom has an energy < THERMAL. 
1 if the site Lis occupied by any atom. 
0 if the site L is unoccupied. 
The number of the time step in which atom L 
exited the lattice. Zero if Lis in the lattice. 
The number of the other atom in an interstitial 
pair with atom L. Zero if Lis not a member 
of an interstitial pair. 
The site occupied by this atom, even if itis 
the original site. If this number is zero, then 
the atom has an energy > ETHRESH and is 
wandering through the lattice. 
The number of the face through which atom L 
exited the lattice, called NFACE. 
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VARIABLE 


OCTAL (L) 


TLAST (L) 


CUROFF™* 


PERCENT* 


USAGE 


OCTAL 
POSITION MEANING 


15 Bit 1 1 if the atom is on the LATER list but has not 
been through INVAC. 
0 if the atom is not on the LATER list. 
Bit 2 1 if the atom has been through INVAC and has 
not had another collision. 
0 if the atom has never been through INVAC, 
or has had a collision since going through 
INVAC. 
Bit 3 Not used 
16 Not used. If Bit 3 is used this results ina 
negative number and should be avoided. 
Equivalent to LB (L). Used for output purposes only. 
The absolute time from the start of the run to the time of 
atom L's last collision. Given in jiffys. 
A percentage of a Time Step Length. All atoms witha 
time of collision within CUTOFF % of a TSL of the mini- 
mum time will take part in a near simultaneous collision 
with atom L. 
A percentage of a Time Step Length. The absolute time 
of completion of an interaction equals (Time at start) 


plus (PERCENT) (TSL). 


on 
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VARIABLE 
TFAC* 


LU LL* 


bg, TY .,..02* 


A* 


SPHI* 


ETHRESH* 


THERMAL* 


FMASS* 


Jeo 


USAGE 
A scaling factor for the Time Step Lengths. 
The interval, in time steps, at which atoms will be shot 
into the lattice, provided the proper jump switches are 
set. Does not affect operation if the jump switches are 
not set. 
The X, Y, and Z fixed point dimensions of the lattice, 
respectively, such as 9x 14x15. IX must be ODD. 
One half the side of a unit cell. For copper this is 
1,807 Angstroms. 
The sphere of influence of an atom. One atom has an 
SPHI of (1.414)(A), and the other atom is considered 
to be a point. Given in Angstroms. 
The displacement threshold energy in electron volts. 
An atom with an energy > ETHRESH may move through 
the lattice. Atoms with energy below ETHRESH are 
forced to remain on their sites in the lattice. 
The thermal energy limit. Atoms with energy above 
THERMAL are allowed to vibrate on their site. Atoms 
with energy below THERMAL are placed on sites, or in 
interstitial pairs by the INVAC subroutine. 
The mass of a lattice atom in Atomic Mass Units (AMU). 
The interval, in time steps, at which the Time Step 
Length is recalculated. 
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VARIABLE 


MNTS* 


MNP* 


FI (1-3) 


AB (1-3) 


NLINE 


NPLANE 


LMAX 


MAXLAT 


TEMP, TEM 


SITE (L) 


TTIME 


OENERGY 


USAGE 
The maximum number of time steps the program will run. 
The maximum number of particles to be shot into the 
lattice. 
Used throughout the program. Usually denotes the 
floating point co-ordinates of a particular site. 
The physical dimensions of the lattice in Angstroms in 
the X, Y, and Z directions, respectively. 
The number of atoms in a line in the X direction for the 
undisturbed lattice. 
The number of atoms in any undisturbed X - Y plane of 
the lattice. 
The maximum, or total, number of atoms. 
The total number of atoms in the lattice, excluding 
those shot into the lattice, i.e. the number of sites. 
Floating point variables used for temporary storage only. 
A short subroutine. Given the number of a site as an 
input, SITE computes the fixed point co-ordinates of the 
site. These are then found in IB (1-3). 
The total elapsed time from the start of the run to the 
end of the current time step, in jiffys. 
The total original energy of the lattice in ev. This is 


increased every time a particle is shot into the lattice. 
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VARIABLE 


ENERGY 


TENERGY 


IENERGY 


MISTAKE (1-10) 


NA 


USAGE 
The energy in ev. of the most energetic atom in the 
lattice. Atoms that have left the lattice are not con- 
sidered in the computation of ENERGY. 
The total energy in ev. in the lattice. This is computed 
every time step for comparison with OENERGY. 
The number of the most energetic atom in the lattice. 
An array of 10 error counters. Their meanings are: 
l The number of time steps for which TENERGY varies 
more than 1% from OENERGY. 
2 The number of times a vacancy is filled, or an 
interstitial is formed, near an edge or corner of 
the lattice. 
3 The number of times a negative energy was computed. 
4 The number of times an atom formed an interstitial 
with — i.e. a false interstitial. 
5 The number of triple interstitials formed. 
6-9 Not used. 
10 Set equal to MISTAKE (2) at the start of INVAC so 
aes the increase in MISTAKE (2) may be calculated. 
The index for the TSLI and NTSC lists. Always one 
greater than the number of times the TSL has been 


recomputed. 
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VARIAB LE USAGE 

NEXT The number of the next particle to be shot into the 
lattice. NEXT is one before any atoms are shot into 
the lattice. 

MINTS The minimum number of time steps, i.e. the lower limit 
of the time step DO loop. Always one unless the re- 
generation method is used. 

T1ll The square of SPHI. 

MOVE An indicator used to determine the exit point from sub- 
routine INVAC. 

NOW The lattice site nearest an atom, computed in INVAC, 
used to determine an atom's oscillatory center, if the 
atom's energy is below ETHRESH. 

LOCATE An indicator used to hold the number in Index Register 5, 
enabling us to use Index Register 5 without destroying 
the contents. 

CON (1-4) Constants used in the subroutine INTER. These have no 
relation to the main program. 

W (1-4) Constants used in subroutine INTER. These have no 
relation to the main program. 

XC (1-6)* Only two used. These are the constant parameters used 
in the Born-Mayer potential, and can be changed to fit 
any one of the three types of Born-Mayer for copper, 
as devised by Gibson, etal. 17 
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VARIABLE 


Cr Gi 
MASK] - 


MASK8 


ALFA 


ALPHA* 


BETA* 


YENTRY (NEXT) 


ZENTRY (NEXT) 


EPB 


TSL 


TSLI (NA) 


USAGE 
Constants used in subroutine INTER. Computed outside 
the subroutine as a time saving measure. 
Masking quantities used for logical arithmetic with 
LB (L). The definition of the significance of the octal 
positions in LB (L) and the masks themselves explain 
their specific applicability. 
107° A small number for comparison to floating point 
variables to prevent use of insignificantly small quanti- 
ties. 
The angle between the original X and Y velocity compo- 
nents of a particle shot into the lattice, given in radians. 
The angle between the original X and Z velocity compo- 
nents of a particle shot into the lattice, given in radians. 
The Y co-ordinate at which a particle shot into the lattice 
crosses the X= 0 plane. Given in Angstroms. 
The Z co-ordinate at which a particle shot into the lattice 
crosses the X = 0 plane. Given in Angstroms. 
The energy in ev. of a particle shot into the lattice. 
The Time Step Length in jiffys. TSL is calculated so that 
the most energetic atom in the lattice will move a dis- 
tance equal to (TFAC) (A) in one time step. 
A list of the first 30 TSL's that are calculated. Used 
for output purposes only. 
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VARIABLE 


NTSC (NA) 


EMAX 


j 


LATER (NZ) 


NZ 


INMIN 


INMAX 


INPUT 


USAGE 
A list, corresponding to the TSLI list, that gives the 
time steps in which the values given in the TSLI table 
were calculated. 
A list of the values of ENERGY for the first 30 time 
steps. Used for output purposes only. 
The sixth Index Register, and the number of the current 
time step. The current time step can always be found 
by stepping the computer console and recording the 


value of this Index Register. N is used for nothing else. 








The number of the atom sequentially under consideration. 
The index of the DO loop starting at Statement 160. 

The number of the atom currently under primary considera- 
tion, usually equal to L. 

Usually the number of the atom considered in relation to M. 
A list of atoms whose energies are less than THERMAL 
and are to be put through the INVAC subroutine. 

The index for the LATER list, used for no other purpose. 
The index for the LAST list, used for no other purpose. 
The maximum number of atoms listed on the LAST list, 
equal to NMAX when the list is made. 

An indicator. If zero, then the KHIT list is copied onto 


the LAST list. If not zero, the LAST list is left untouched. 


97 





VARIABLE 


NUM ( ) 


KHIT ( ) 


NMAX 


TMIN ( ) 


DSTANCE ( ) 


TeX, °) 


ITEMP, JTEMP, 
MTEMP, NTEMP 


NPAIR 


USAGE 
A list of all the atoms in a box of side length 3.02 (A) 
centered at atom M. Only atoms on this list are con- 
sidered as possible targets for M. M is first on the 
list. 
The total number of atoms listed on the NUM list. 
A list of atoms, derived from the NUM list, that atom M 
will hit. Atom M is first on the list. 
The total number of atoms listed on the KHIT list. 
The time, in jiffys, of closest approach between atom M 
and the corresponding atom on the NUM list. This 
assumes no deviation from straight paths. Measured 
with the start of the current time step as a zero reference. 
The distance, in Angstroms, between atoms when at their 
point of closest approach as determined by their TMIN. 
The time, in jiffys, at which the two atoms are SPHI 
distance apart. Will be less than the corresponding 
TMIN if there is to be a collision between M and the 
atom on the corresponding NUM list. 
Fixed point variables used throughout the program for 
temporary storage purposes. No particular significance. 
An indicator for a computed GO TO statement. NPAIR is 
3 if one of the atoms being hit is a member of an inter- 
stitial pair. 
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VARIAB LE 


IN3 


LAST ( ) 


Tl - T10 


NFACE 


JA 


KMIN 


NHIT 


ELOST 


USAGE 
An indicator that shows the number of those on the NUM 
list that will be hit by atom M. 
A duplicate of the first KHIT list comprised for M when 
M=L. All atoms on the list will become M's when the 
preceeding M is no longer able to have collisions in this 
time step. 
Temporary variables (see BOX 17, APPENDIX IV) used to 
speed up computation time. 
The face of exit of an atom when it leaves the lattice. 
Faces 1, 3, 5 are the negative X, Y, and Z faces, res- 
pectively, while faces 2, 4, 6 are the positive X, Y, 
and Z faces respectively. 
The index for the KHIT list that shows which atom 
atom M is colliding with at the present time. JA's range 
is from 2 to NMAX. Has other uses in INVAC. 
The number on the T list that is smallest. T (KMIN) is 
the minimum time, and is the collision time for atom 
NUM (KMIN). 
The number of atoms on the KHIT list that have been 
hit by M. 
The energy in ev. that M lost in its collisions with the 


other atoms on the KHIT list. 
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VARIABLE 


SAVE (I,K) 


GOLD ( ) 


PAR 


TIME 


DEV (I, 1-2) 


USAGE 
The pre-collision velocity components and energies of 
the atoms on the KHIT list. SAVE (1,1) is the initial 
X velocity of KHIT (1). SAVE (2,1) is the initial Y 
velocity of KHIT (1), etc. SAVE (4,1) is the energy of 
KHIT (1). 
A floating point variable used in the energy and velocity 
scaling process. COLD (1-6) correspond to the first 
six co-ordinates in the B array. They are the sums of 
the individual components of the DEV array after each 
collision, i.e. COLD (1) is the sum of all the DEV 
(1,1)'s for a series of collisions. COLD (1) is then 
the sum of all the X co-ordinates M would have had as 
a result of the collisions with the rest of the KHIT list. 
The scattering parameters in meters, and the third input 
to the INTER subroutine. 
The time T for M's previous collision with one of the 
others on the KHIT list. Forinstance, if JAis 3, we 
are considering KHIT (3) in relation to M. TIME is then 
the T corresponding to KHIT (2). 
An array, used in subroutine INTER, holding the position 
co-ordinates, velocities, and energies of the two atoms 
at the end of the interaction. The first index, I, has 
the same meaning as the first index of the B array. 
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VARIAB LE 


STUPID 
SCALE ( ) 
TURN* 

IT 

MANY, NANY 
JPMAX 
MOON ( ) 
NOW1, NOW 


USAGE 
The second index is one for atom M, and two for atom J. 
The B array is not changed in the subroutine, but the 
DEV array contains the output of the subroutine. The 
B array is then scaled as necessary. 
A scale factor used in the velocity scaling process. 
A three member array used in the velocity scaling process. 
An atom in the intermediate energy range must lose TURN % 
of its energy before it is re-directed back toward its 
oscillatory center. 
A temporary variable used to facilitate programming. Also 
the index for the DO loop over the LATER list. 
The atoms occupying the site chosen for atom M after M's 
energy has fallen below ETHRESH. If NANY is zero, then 
M and MANY will become an interstitial pair. If both are 
non-zero, a "triple interstitial" is formed, while if both 
are zero, M occupies the site alone. 
The number of atoms on the MOON list, i.e. the number of 
nearest neighbor sites. 
A list of nearest neighbor sites for a particular atom. 
Temporary variables used to indicate the site chosen as 
oscillation center for an atom in the intermediate energy 
range. NOW] = (NOW shifted 27 binary positions to 
the left). 
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VARIABLE 


NAVY, MAD 


ZZ (6) 


NN, IL, LL 


MJ 


MM 


IFACE ( ) 


NOON ( ) 


USAGE 
Similar to LOCATE, except these are used for Index 
Registers 2 & 1, respectively. 
The time in seconds necessary for the interaction. 
Variables used to denote atom numbers. NN is the 
index for the DO loop starting at Statement 241. IL 
usually denotes a site number, while LL usually de- 
notes the number of the "other" atom in an interstitial 
pair. 
The number of the atom used as an input to INVAC. 
Usually equals MJ, but if MJ finds no vacancies then 
MM is the prospective interstitial partner atom. 
A three member array denoting the face an atom entering 
INVAC is near (i.e. within one lattice unit). The index 
is 1, 2, or 3 for the X, Y, or Z directions, respectively. 
The list of vacancies near atom MM (which is usually MJ) . 
This list is derived by examination of the MOON list, 
i.e. if an atom is not on the MOON list, it will never 
be on the NOON list. 
Either 1, 2, or 3. Denotes the direction (X, Y, or Z re- 
spectively) of the axis of the prospective interstitial 


pair. 


NOTE: This is the end of the list of variables for the MASTER program. 
Many of the temporary indices are not listed in this table. 
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The symbols listed below are peculiar to program SLAVE and the 


usage given applies only to program SLAVE: 


VARIABLE 
NTS 


KET (1-26) 


IHISTO (1-105) 


NBUL 


EBUL 
NSUBST 
NSAME 
NVAC 
INTERST 


ETHERM 


NSIDE (1-6) 


ESIDE (1-6) 


ECAL 


LMAXCAL 


YAN GLE 


USAGE 
The actual number of time steps the run completed. 
A list of the energy ranges for the histograph distribu- 
tions. 
The histograph format array. 
The number of knock-ons (i.e. atoms with energy above 
ETHRESH) in the lattice. 
The total energy of all knock-ons in the lattice. 
The number of replacements in the lattice. 
The number of atoms on their original lattice sites. 
The number of vacant lattice sites. 
The number of interstitial pairs in the lattice. 
The total energy of the atoms that have gone through 
INVAC. 
The total number of atoms that have left the lattice 
through each of the six faces. 
The total energy of the atoms that have left the lattice 
through each of the six faces. 
The same as TENERGY, except calculated by Program 
SLAVE. 
The same as LMAX, except calculated by Program SLAVE. 
Angle ALPHA given in degrees rather than in radians. 
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VARIABLE 


ZAN GLE 


USAGE 


Angle BETA given in degrees rather than in radians. 


NUMBERS (1-50) A fixed point array used for temporary storage only. 


IC (1-100) 


DEPTH 


PERDIST 


EL, 
XYZ (1-3) 
NTOTAL 
NU MBIG 
TOTALN 


BIGNUM 


A fixed point array used for temporary storage of the 
histograph channels. 

The total penetration of an atom in the positive X 
direction. 

The total radial penetration of an atom perpendicular 
to the X direction. 

Usually the number of an atom, temporary usage only. 
Corresponds to the FI array in program MASTER. 

The sum of all the numbers on the NUMBERS list. 
The largest number on the NUMBERS lisi. 

NTOTAL in floating point format. 


NUMBIG in floating point format. 


The symbols listed below are peculiar to program RON and the usage 


given applies only to program RON: 


VARIABLE 


T 


DX 


- USAGE 

Temporary floating point variable. Not to be confused 
with the T array in program MASTER. 

A distance parameter, in Angstroms. If any of an atom's 
co-ordinates differ from its stored values by DX, the 
atom's position will be printed, and the new position 


will be stored. 
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VARIAB LE 


P (1-3) 


IN 


X (1-3,N) 


NT (I) 


NO 


S (1-3) 


USAGE 
Reciprocal of A. 
The co-ordinates of the atom under consideration. 
The number of the time step the atom under considera- 
tion has co-ordinates P (1-3). 
Stored co-ordinates of an atom to be written out at a 
later time. 
An array used in the output section. The lower half of 
a storage cell is the number of the atom whose co- 
ordinates are X (1-3,N), and the upper half is the time 
step the atom had these co-ordinates. 
An indicator. If greater than zero, a heading has been 
written on the output tape for the atom under considera- 
tion. If less than zero, no heading has been written 
on the output tape for the atom under consideration. 
Co-ordinates of an atom in floating point form in A units. 


One A unit is 1.807 Angstroms for copper. (See A, above) 
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LATTICE CONSTRUCTION 
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Z =(77-1)/NPLANE = 76/120 = 0 Remainder = 


Y =(76) /NLINE 
X =(4)(2) = 8 
+Z=9+0 


= 76/8 = 9 


which is ODD, therefore X = +12=9 
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APPENDIX IV 
PROGRAM FLOW CHARTS 


SYMBOLS USED IN THIS APPENDIX 


START of a program , subroutine , box , er 
Pegit 


Read in or write out 


Decision or branch 


General other computer operation 


Position desiqnation 


Logical product of IJK and LMN 
Contents ot the “A register \ 


Contents of the "“Q register 


Statement number 132 
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PROGRAM MASTER Boxes I|-5 















READ and WRITE 


Pive comment cards 
(tape 2 to tape 3) 


WRITE: “THIS 
TS PROGRAM 
MASTER™ 








READ Ctape 2): CUTOFF, PERCENT, 






ae TEAC) LULL LY, IZ, SPAT, ETHRESH, 
- THERMAL, FMASS, JPB, MNTS, MNP 
NLINE= 
FIQ\)=Ix 
Rut 2)= BY | YMOM(I-3)=O | START {| Ex+/2 
mis)- lz AB(\-3)=A* FI(I-3) Box 4 NPLANE = 





NLINE «(IY+1) 





DO 1020 START MAXLAT = NPLANE -CiIz+1) 
L=!,LMAX Box 5 LMAX:= MAXLAT 


WRITE (tape ony 

TX, TY, TZ, AB0=3), 

MAXLAT, A,SPAI , FMASS, 
— ETHRESH, THERMAL , MNP 


Compute X,Y and Z 
coordinates of each site 
by usin subroutine SITE(L) 
Then puts L in (10-13) octal 
positions of LBCL), puts 

a “Lin Sed bit of LB(L), WRITE Clape CG): 

(BLI-3, 5) for T= 1, MAXLAT), 
NUINE, NPLANE, YENTRY, ZENTRY 


and sets TLAST(L)=-10.0 


DO 1020 T=4.6 
Set the velocity 1B] 
equal to zero, Set energy 


equal to zero, 
Box 
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Boxes 6-10 


RT Set TTIME, OENERGY, ENERGY, Z2(4), MISTAKE, START 
ae all equal to zero, Set NEXT, NA oe 
Box 6 MINTS equal to one. Ti\=(SPHI ) Box / 





Set MASK’'s for logical Set initial 
START 
ith LBCL), and constants for 
——— Box 8 interaction subyoutine 





explanation of sense switches 





At this point the veal problem starts, all revious to 





a! this was in setting up the problem. The flow 
a chart will be move detailed from here on out. 
READ (tape 7): All ue START 
Reaeneration \nput Box \0 

DO %000 
EWIN 
eS ms Mt dat 
Box 


\1 











WRITE (tapes 344). 
TENERGY, NN, TENERGN, 
ENERGY 











1+ 3 down 


Z00| WRITE (tape %): Yes 
output foe SLAVE 
2 up 



































ia 
Switches 
\43 up 
2 down 





WRITE (tape 3): 
BCD part of 
output 


Stop it 
Slee daltalh 
1is set 









Stop if 
Stop switch 


\ ig set 





Com pute 
(N+I)7 LULL 


Go To 4999 
(Abandon the run. ) 


ZO GoTo ~ ZO 
(20* 









FO 










ENERGY- 
THERMAL 





>O a0 
GoTo 
800) 
we el 
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LMAX =LMAX +1 

Read in AY,2, Samoa 
angles x, B 

Vy = Vy tan of 

Ve = Nig tan B 
LB(LMAX)= 2 
TLASTCLMAXN) = - 2.0 
Y= "4 - V Ny 
Z=Z2- Vy 


TSLICNAzTSL 
NTSC(NA)= N 
NA =NA+ | 





If TENERGY 


1s within vanae 


TYOENERGY — |.0\ OENERGY 






Boxes \2-15 











OENERGY = OENERGY + BC7,LMAX) 
TENERGY = TENERGY+ BC7,LMAX) 
YENTRYCNEKT) = B62, LEMAX) 
ZENTRYCNEXKT) = 8 C3, LMAX) 
NEXT =NEXT +1 

EPB= (7, LMAX) 

Recalculate TSL, store in 
TSLICNA), Store N in NTSCCNA) 


START Compute 
Box (2 N/JPB 
Recalculate TSL = 
on basis of ENERGY 
+0 
TTIME = START 
TTIME+TSL | “\ Box \3 


START Add one to DO 160 Add 8(7,L) to TENERGY 
Box \4 MISTAKE (\) L=1,LMAX Clear MASKS bits of LB(L) 





> 


EMAX(N) = 
ENERGY 





TENERGY =o CO 
ENERGN =O 


Rox 
Sy 
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Boxes 16-]/ 


Go To 
$000 


#O 









Et 


TLB(L)-777716) <amd> 

NUM (1-20) = O 
KHITC-20)=0 ==) 
LAST (i-20) =o 
T(1-20)=0 
DSTANCE (1-20) =o Ca 
TMIN(I-20) = 0 MeL L a INZ=0 

Box I6 





J-LMAX | LB(3) TIT718} 





“<B> 
30 

Go To 277 

5 (Box 25) 


NUMC(K)= 3 
K=K+I 





K=K-I a 20 
KMAX=K ee 


<0 


Go To 277 
(Box 25) 
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Box 17 
DO 215 ITEMP= NUM(K) 
K= 1, KMAX 
TT = T474TS7+T6* LL 
Te = Vi2+72°+T3" 








Ti= BC), ITEMP) - BCI, M) 
T2= BC2,ITEMP) - 8(2,M) 
T3=B(3,ITEMP) — B(3,M) 
T4= BG ITEMP) -BG,M) 
TS= BCS, TTEMP) -B(S,M) 
T6= 8(6, ITEM?) -BLE,M) 







Oe liGis Te 
D4 S65 Z 
OSTANCE(K) - 
20 SPHI 


T9=T1-T44T2:T5 <0 
eS aC OSTANCE(K) = {LTI+T4 -TMINGRY* 
TMIN(K)= - t + FT 2+ TSATMINCKY] 7+ (£3476: TMINCK)]~ 
{ 
214 | Le 
ZO 
CC tae SL 
Kox | <0 
\g 

>O if Compute T(x) +TTIME. 
CA) ~TSL— TLASTCITEMP) 


=O 
> 


TK) + TTIME A 
Ope} INZ=INZ + | 2Q\5 
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TK\= TMIN(K) 













eo 
T(R)+TSL 





—TSL-TLASTCM) 








Box \8 
#0 JA = KMAX ~| 
TEMP= (000: TSL 
apa) 


=O 
GoTo 277 DO 217 
(Box 25) K=2, KMAX 
K =K+l 


L=KMIN 
oi 2 |] K=O =a 
Tl=CUTOFF-TSL 


“9 w 
Be a el 
K=|.KMAX TCK)-T(X)-TI <a KMIN =K 
’ 


| 
Go To | KHITCIT) =NUM CR) y=T-\ 
J=I9+\ NMAX= 5 


INMAX = J 
VO 300 ye 
LAST(K)=KHIT(K) Met 
=, J 
#0 


EC | INPUT = 
Go Vo 277 (Box 25 <reed> — 30 
#0 
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Bo x | 7 


J =KWITCK) 
O 35) 
2 SAVE (1-3, K)= COLD(I-7)=0 
K= |,NMAX B(4-6, 3) 


Tf 
Determine K so that =O [Le ecs)-MAsk J =KHI(JA) 
NUMCR) is J hee : 


PAR= DSTANCE(K): 1907/2 (A ) 
NPAITR= 3 





TIME =TCK) 
CALL INTER(M,S, PAR) 
TLASTCS) =TTIME -TSL +TIME 





Move #J ahead to time TK). 
If JA=2, move *M to TCK). 
Tf TA#2, move 44M back 

From last T(K),and then move 
ahead to time of present TK) 






Add DEVA-6,1) to 
COLD (1-6) 


(1-7, T= DEV(I=7,2) 


(A ) NHIT=NHIT+1 | Add £8(7,3) 





- SAVE(4, A} to ELOST 


=0 
Box 
20 


no? 





Go To 
Box 2| 
<O 
> 


STUPID = 


VTEMP/TEM 


A) Biéx 20 


TEMP = 0 
TEM=0 


B(4-7, M) 
=O 









TEM=I/NUHIT 
BC(\-3,M) = 
TEM - COLD(I-3) 










y= KHITCK) 
B(4-¢,3)= 
STUPID BG-,7) 
Saino) = 
STUPID: B(7,5) 













J =KHIT(K) 
‘Add = SAVE(I-3,KY* 
to TEMP 

Add = B(4-6,5)~ 
to TEM 






TEMP= S CoLD(4-6)~ 





2 snes iies 
feo 6518 FMASS TEMP 


B(4-6,M)=TEMP> COLD(4-6) 









B(4-6,M)= 
DEV(4-6,1) 


CALL SITECL) 

SCALE (I~3) = 
A-TR(I-3)=B(I-3,m) | 

TEM= = SCALE(I-3)? 








T URN-SAVEG,1) 








T= [Le()-Maska}-y 
TEM=0 










TEM= T2/G5IQ«+ FMASS TEM) 
B(4-6,.M)= TEM-SCALE (1-3) 











TLAST(M) = TTIME-TSL4TIME 
+22(6)-10'7. PERCENT 

TEMP = B(7,M) 

LOCATE =O 







T= [Le(m)-magKd] «2° 
Cleay 3rd bit 
of LB(I) 


Go To 513| _ #0 
(another page) Clear bit #2 of LBM) 
Clear MASK2 


LTEMP= 10-13 octal positions 
TTEMP= 6-4 octal positions bits of LB (M) 


TTEMP + Legg FO =0 =O) Si Ze Set bit #2 
YTEMP>2 of LBUM) 


Clear MASK3 bits 


If 


=0 
(A) LB(M)-MASK2 










Clear bit #3 of Rent eee Cleary MASKI bits Go To S\4 
ITEMP= 0: o 
of LB(ITE MP) LBCM). Clear of LB (ITEMP) (another page 





bit 43 of LBM) 


7 
LATER(NZ)= M a = se 
NZ=NZ+ | pis 
Page 
Clear MASK7 
MISTAKE(3) = bits of 
MISTAKE(C3) + | LBCM) 


MJ =M 
LOCATE =O 600 MOVE = | 


NOW =O 
CALL TWAC (MT) 
MOVE =O 
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th 
m 








FO 
=O 
= Go To SOT. 









If 
fLe(r)- MASKI | 


-NOW| 






*© Next page 





DO 507 
I> 5 LMAX 


’ : 

NANY = Cleay bit 2 
\, 

Add one to =O 

MISTAKE(S) = 


Put ITEMP in 
LB(m) 

Put JITEMP in 

LB(NAWY) +LBCr) 


Now = Put NOW! in octal 
Now-277 positions !0-13 of LB(M) 


ITEMP=I-2"° 
JTEMP=M-2!° 


mEO 


=O 


LBCNOW) 


Go To 513 
(next page) 





DO Sil 
K=|,JPMAX 


< 
— eT 
DO 512 
DSTANCE (kK) 


=K 
>O —~ DSTANCE(Z) 


NOW = NOON(T) 





Now \ = 


Now: 2!> 





Go To 
Box 22 


of 
preceeding page 


IL= NOON (Kk) 
CALL SITECILY eal 
Add [RC1-3,M)-A-Ie(1-3)] to DSTANCE 








K=2,TPMAX teu 








(A) T2= BM) 


pase 
22 BS a 


J) = NAVY 


Add one to 


M= LOCATE 
T=MAD 


Go To 235 
(Box 22) 


MISTAKE (3) 








/ nat 
(sw ee 


Clear MASK2 bits of ITEMP = octal \0~13 >oO 
LBCJ). Set 2nd positions of LB(J) 
bit of LBCI) in upper half 








ae 
Set bit 2 of =O 


L BCT) 






Go To [L@(3): MASK} 
233 
Ae) 
Clear 3rd bit 


of LBC ITEmP) 


Go To €00 
(Box 21) 
LATER (NZ) Cleary MASKS 
=f and set MASKG 
NZ=N2+| of LBCS) 


b 
A 


Clear bit 2 of LACD 


LDA2(Le), LDQ(MASK2), 
STL CITEmP), LDQ(MASK)) 
STL (ITTEMPY 









LOA(XTEMP), ARS) 


STA(ITEMP) 
K= I TEMP 
Clear 4th bit 6 LAR) 


Clear MASKI bits ARS(15) 
of LBC J TEMP) STA CJ TEMP) 












Caina> =O Box 
me 
#0 
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f, 


= Boxes 23-25 
START =O o To 
cas 


70 


DO 2€2 DO 20 
J=KHI T(NN) 
NN =I, NMAX K=1,KMAX 
B(\-3, 5) = B1I-3, 3) 


£0 
- [TIME +2206):10° J+ B(4-6, 5) 





De te Box 25 Lor K=2,NMAX , D=KHIT(KS 


START =2 | From LACT) find if T was part of an 
\nterctitial pair. If co, put J's partner 
Box 24 - on the site of the interctitial pair . Clear 










#3 the site number and partner number Leom 
LBCJ). Also clear the partner number from 
START Leary), If JT was raf part of an 
277 Bax as interstitial pair, a0 to Box 255 


Go To 200 


Tf 
LB(M):MAS 
Box Ib [LB(M)-MASKG} 


INMIN = 
INMIN + \ 


Lik3 





Boxes 26-28 


START meer 
K =LAST()) 
Box 6 Kousreo BCI-3, K) r BI-3, K) 
+8(4-6K)> TSL 
Set bit 4) and clear 
MASKE of LBCK) 


My = LATERCIT) CALL INVAC CMJ) 
Upper half of IU =N 
Lowey hal€ &F TU =MIT 


Write Tape 6 , IU, & (i-3,M7) 


: 





rf 


ENERGY - 
BC7,K) 






ENERGY= BC 7,K)| . <9 
TENERGY= Kk. 


















START 
Box 27 





GE 










79 E DO 261 a 
If 
START DO 7998 DO 206 fm 













DO 252 
=. Ss 












B(L,M)-AB(I) 


Go To 
M= LAST 
St) 1998 








NFACE=2-T 


NFACE=2-I-| 





ce 


S000 





Put NEACE in octal #14 of LRLM) 





ITEMP= octal [0-13 of LBCM) 


Box 
1998 NFACE = 0 ITEMP= octal 6-9 of LB(M) 
a 7498 weace=0 Cleary MASK3 bits of LB(M) 
ds Clear 48 bit 4#0_ | Clear MASK\ bits Qo 
of LRCLTEMP) of LB(JTEMP) 
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O 





Box 29 


ee DO 9000 #0 
STUG(IT) L=1, LMAX 


e 0 
Bree of LBCL) BC\-3, b) 


“293 
. O - BiI-3,L) = B¢ ) 
= =O I-3,L)= 60-304 
[LB(L)- 777718} (LB (L)-MASK7] 
+8(4-6,L):TSL 
FO 
Go To Go To BIEL = DO 312 
4000 320 AB T=1,3 
| ae >O : 
Go To 
312 NFACE=2-J 
een> — [wacesa 
<O 


NFACE=2:I-| 






Store NFACE 


and N in ITEMC= octal 10-13 : 
B(L) JTEMP=octal 6-4 Clear MASKS bits of 
of LR(L) | BCL) 






Cleae 4B =O 
LBCLTEMP) 


Clear MASK| bits of LR(TTEME) 


——_ 









<O 






re 
ENERGY— 
BL) 


4000 ENERGY = BC7,L) (B) 
Box SO TENERGY=L : 
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Box 30 






WRITE tape 6: IX, IY, IZ, ABCI-3), NLINE, NPLANE, 
MAXLAT, LMAX, A, SPHI, FMASS, ETHRESH, THERMAL, EPB, 
LULL , ALPHA, BETA,N,MNTS, TTIME, TSLI NTSC, 


OENERGY, TENERGY, MISTAKE, YENTRY, ZENTRY, NA, 
EMAX, EFOUND. WRITE tape €: CLBCL), B(I-7,L) L=1,LMAX 











Go To 
49498 









Tf 
Switches 


2. Down 
\,3up 


Yes 





Stop if Stop 
Switch | is up 






REWIND 8 






No 


WRITE tape 3: TTIME , TSL, ENERGY, N, (or MNTS), TENERGY, 
OENERGY, MISTAKEC1-5) , TX, LY, LTZ,b , BC1-7,L) , OC TAL CL) 
GoTo II6 Stop it Stop £0 

499 
(Box I) Switch lis up 


=6 


WRITE tape 3:3 


“THE RUN WENT THE 
FULL DISTANCE * 









WRITE tape 6. 
oo 





END 


REWIND 6 }- 
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SUBROUTINE SITEC(M) Box 3] 


cag ikea 
(0 zt TBS) = NeLANE 
\ 
IB(=vemainder * 2 eae tannin ons 


Odd 









T8(2)+IB(3) IB(N= IBCs | 


END 


Pee 





SUBROUTINE INVACCMJ) Box 32 


Cleary MASK2 4MASKG , set 
MASK7T of LBCMM) 
















NIF =O 
Go To DO 2370 

IFACEL)=0 3 MISTAKE (io)= 

eee Be MISTAKE(2) 
20 

>O a £O | IFACECTD= 2:2 Go Yo 

B(z,MM) NIF =NIF+! 22370 
<0 


IFACE(I)= 2-I -1 YO 2337 TB(1) = B(IT,MM)/A 
TB(1) = IB(r) +1 










B(T,MM) 7A 
-TEMP - 0.5 





Box 
so 
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Box 33 







Even 


TL= 183): NPLANE + 
TB(2). NLINE+I8(0/2 +\ TGC!) +IBQ2)+IB(3) JM= | 
FO JM= 
=O 


TL= O beh Odd ¢ 
TA=0 <r LO. I 


TRE)+IB(D 
SL=-| 


DO 2115 | Tf >O JAz= JA+tl 2B 


JR= TFACE(L) (BY 
3K 


= omg ED 


(A) >| MISTAKE (2) = MISTAKE(2) +1 ye 


IM:6'T0 
+1+ 78 
JG 


<4 +2 1S . 
SPMAX= 5 <Gw> <2 
ES 





















NOON(I-10)=0 
MOON(I-\5)=0 
t= \ 






or 





L1lg 





Box 34 


DO 2475 Go to statement 


JP=|,SPMAX number (2700+3W) 


270 | GO TO €2455,-2451, 2452 ,2453, 2454, 2445), TP 





2702 GO TO (2455, 2451 , 2470 ,2456, 2452, 2453, 2457, 
2458 . 2459, 2454, 2460, 2461 .2462),5P 


2703 GO TO (2455, 2451 , 2452,2453, 2454), 5P 





KO 
at 
© 
3 


GO TO (2451, 2452, 2453,2454, 2465), IP 


105 GO TO (2455, 2451, 2453, 2454,2465), JP 


Wo fo 
J 
© 
oN 


GO TO (2455,2452, 2453, 2454,2465),IP 





2107 GO TO (2455, 2451, 2452, 2453 ,2465),5P 


210% GO TO (2455, 2451,2452, 2454, 2465), IP 


2704 GO TO (2455, 2410, 2456, 2457, 2458, 2459, 2460, 
2461,2462),7P 


2710 GO TO (2455,2451, 2452, 2453, 2454, 2459. 2459, 
2461 ,2462).,7P 


Zip GO TO (24556, 2451, 2470, 245932, 2457, 24538, 2454, 
2460, 2461), 7P 


GO TO (2455, 2452, 2453,2454,2456 , 2457, 2454, 


22 | 
IAeS, 24:40) , SP 


GO TO (24565, 24-51, 2452, 2453, 2456 , 2457-2458, 





2459, 2470 D,TP 





GO -TO -(2.4555-2451.,.2470, 2456, 2452,.2454.-2460, 
2461, 2462), IP 





120 





=O 
(A) 






CA) will contain 
245) XL+NLINE+TL. 
2452 IL=NLINE +TL 
2453 IL + NPLANE?TL 
2454 IL- NPLANE +JL__ 
2455 TL ei, ey 
2456 JL~NLINE+IL+L 
2457 IL+NPLANE+SL41 
2456 IL+NPLANE*+NLINE 
2459 IL+NPLANE-NLINE 
2460 IL -NPLANE+JL+1. 
246\ IL- NPLANE *+NLINE 
2462 TL- NPLANE ~NLINE 
2465 fe) — | ee 
2410 IL+NLINE+JL + | 


) START \ 
Go To 274\| 40 =O DO 2485 
runt WEG 


NOONCTP)= L. rf r= 
Bex FT. 
36 






L# 







MAXLAT 






<O 






Moon(r)= ITTEMP 
Tg | 



















hn 





r <0 . 
START 3 Boxes Sb+3/ 
Box 26 =O ZO 

>0 =O 


IL= NOONG) GoTo 27 
BTS ea gy] LELENQONG) | [Go To 2000 oT> 2720 


@ALL <2s7EGO) |( eox888) (Box 3%) 
2442 


= i IL) tnto 
B(1-3,MM)= TBCI-3)°A NI=3 Put new site ATL) tnt 
LBCMM), and set 4B in LBCL) 


RETURN 
VOWS 12 
= | WRITE (output). “ATOM =O 
os a | CMM) has es placed on 
Site sai ) 
#0 
LTEMP = NOON(K) DSTANCE(K) = Go lo 2800 
DSTANCE(K)=O (ene 23) 
CALL SITE (ITEM?) z [TB(1-3)-A- BC-3,mi) 
J=i DO 25 (14 
K=1,7* 


M=\| 
MOON (M) = 
NOON (3) 







DSTAN CE (rk) 
~ DSTANCE(3) 





MOOoN(K)= 


£ 
DO Adib ee 
Kel yoK 
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ee. =a Ts 
Gm) ED <B> [Be 


Go To 2492 
(Box 36) 






IL= 
MooNn(K) 











yu = 
MOON(2) 






TL=MOON Z\ ) 


T1=Z BG-,mMy" ITEMP = MOON(R) 
TA | DO 2645 TTEMP= = 
K= |, SPMAX TTEMP. 22 
TMIN(K)= | _ #0 
1000 'TSL 


=O 


Go To CALL 
#O 2 = Tae+ 
ET (i-3)= [8C1-3,mm-ert-d] *BC4-6,MM) 
A-TR(I-3) 


TMINCK)= 
Add [BCI-3,MM)-FIC- to TMIN(R) -T2/T\ 





DO 1700 
T=|,LMAX 










Go To KMIN =| 
Next Page TK = TPMAX - | 
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neCe ) = 
A&SE ( 
B(I4MM)-A-I¥(2) 


GoTo 236I 
(Box 32) 





P= 


Moon) CK MIN) 
CALL SITE (re) 


266S 


DO 2680 
T= 1,LMAX 


Does 


T ocey 
sitet ‘4 














Too 
MAX= | 





< 
se 


MAX 


KMIN 


T=(3JB+1) +2 


=K+\I 
pO 26715 
T=l,2 


y 


=T+| 


TE=rY 
1IN2=0 


Go To 2492 
(Bex 36) 


Go To 
Box 34 





Go To 2492 Box 54 


(Box 36) 


TEMP = CALL DO 2725 
BMAX, MM) — BlmAX, MS.) eaneen 














TZ = PSG) 
BCI,MM) =A+T2 
B(IT,M3)= A-T2 






Subtract 4A 
Lom B(MAX, MM), 
Add 4A to 

B (MAX, MT) 





<O 
<> 


Add 6A to BMAX mm) 
Subtract AA fron BCMAK,MT) 










Store LP (site) in LBCMM) and LB(MS) 


Set (4B) tn LBCIP) 
Store MT in MASK bits of LBC(MM) 
Store MM in MASKI bits of LR(MT) 


Z_ Add ene to MISTAKE (4) 









SUBROUTINE INTER(MQ,II,S) 


START XP(\-6) = BCl-6,MQ)-B(I-€,T2) 
Box 40 XCCI-6) = Ve XPCI-6) 


Z20) = C2/exPe(zx7)¢3 
22(8) = Z2(10)- [22(2) -C€9] 


Z2(2) 


Box 41 | 4 | 


= C3+22( 2) 22(id)> L227)" 
22(05)= 22()- Ci -220¢8)] 


Go To 
1S 


20 
22(1)=22{2) 
meet) = 22(2)~ 


ZzZz(%)= C2/ Expr [e2(D-c3] 


224) = 


\— 22(2)/C\ 









2 2(2)+|.000001 
-22(7) 


Z2(5)=O 
Z22(6)=0 
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I? 


Boxeg “+0-41 





Z4(4)=(Relative velocity) ™ 
Z7(10)= 2/CRelative energy) 
Z2(IV= Relative velocity 
z201)= $ 

22(\3)= Anole between line 
jotning the two atoms and 
relative velocity 

22¢(7)= Cl Cinitial quess 


at closest distance of 
a 















2z02) = Z22(C7) 
| - 22(2)- 2205) - 2200) 
2.0: 22(5)+22Q) 


222) is the next quess 
Gor closest distance of approach 









<0 


720) = 
Z2C7) 


Next 
Paae 





Boxes 42-8 





22(9) = 22(4)* CON(IY 


72(2)= Cz2(ayy7 
2201S )= 22(3)+2200/274) 
22(\4)= Z2CN)-£ 2-224] 


ZZ) = CI /ZEXPE[Z20)-22() -¢3|-C9 7Zz012)= \ /[ 1-224] 


Z2(4) = W(T)/SaeTr { ABS L 2209+ 2205): 2208) -2209) | i 





tn | DO 25 
mue T=1,4 






Z2(5) = 2205) + 2244) 
22(6)= Z2(6)+ 224) [2e2ciay] - 


22 (7)=C6*SINELE2(5)] 


; OTART 22()= 2:SQeTFllzz@)| 
hae, Coe eat \ Box 42 72(5)= S* 2X5) 2207) 
3) = 22(5)+ 2213) 22(6)= 22( 3):22(6):207)/220)) 













72(4) = FICI-3N= KeCH-39 X XPLA-O) 





SINE [2%C3){- ZZ())° 3.5 ao |} #2(1-3) = FIsi-35 x XPCI— 3S 
ZB(IO) = 220-3 = 220-38/|200-5| 






CS* 2201) ° COSFL2z2(3)] vector products 





START 
Box 44 Alex t 
Page 


D7 





| Box 45 
O 
on 85) "| eens 
=\, Z22(5)= XC(I)-2ZZ2(8)+ 2Z22(1)-22(7) 


Z2(9) = BCILIT) + XC) + | (23, II) +xXC(I3y]+ 226) + 10 e 


DEV(r,\)= 2209)+ 2206) ZZC5) = ¥C(I): 22.010) + 22CL)-Z 24) 
DEV(I,2)= 7Zz(4)~ 2205) ZZ (9) = BCL3, ITY) + xC (I3) 


71-0 ) DEV(I3Z3,\)= 2209) +2205) 
ER DEV(I3Z,2)= 22094) =2205) 





D0 29 Tienmle £pev(r ie 
Mm 24 





DEV(7,))= FMASS- Tie 0.S51800761F 
DEVC7T,2)= FMASS+T2:0.518001615 








PROGRAM SLAVE Box | 


(sr REWIND TAPE 8 







READ TAPE 8 
TX, TY,TZ, ABC(I- 3), MLINE, NPLANE, MAXLAT ,LMAX, 
A, SPHI, FMASS, ETHRESH, THERMAL, EPB , LULL, 
ALPHA, BETA , NTS MNTS , TTIME , TSLIC)- —“26E 
hese (|20))), OENERGY , TENERGY , MISTAKE (]-10), 
YENTRY (1=10), ZENTRY (1-10), NA , EMAXC 1-30), 
EFOUND, LBC I-LMAX), BCI- 7 \- LMAX) 


Tupe out reminder to 
sel ccldtive switch 2 up 
\ 
READ | CARD (TAPE 2) a eg 
~ U4 


PRINT 4 LINE ) 400 











(COMMENT READ IN) 


NEXT Fi. | 
= KET(26)=1000F READ LN KETC 1-25) 








PRINT 





If 
Switch UP “SELECTIVE JUMP 
2 SWITCH 2 WAS uP 





THISTOCN) = SHCL3,F 
eee a7... 

THISTOC3)= SH3H P. 
THISTO(C4 = 6H4HC. 
THISTOCIOS)=\1H) 









T2QE15 = 327683 


I2ea0= 125s 

T2624 = TQS x Si2 

P2E24= TQE24 

MASK25= 77770B 

MASK64 = 7777000008 

MASK13 = 1777 0000000008 

MASK I4= 7 OCOOOO0OO0COC00COCOB 
MASKSA = 10 COCOCOOOCO0OOR 
MASKS U = 2x MASKSA 

MASKSE T= 73000000C000000078 


-*, 
2) 















NISUBST=0 
EBUL=0O 


NBUL = 0 


: 








CLEAR LBCL) for 
bits in. MASK 64 
and MASK\3 





| VTS=MNTS 





ETHERM = EFOUND 
ESIDE C1-6)=0 
NSIDECI-6) = O 


INTERST=0 
NAA= NA - | 


TSLI(T)= TSLT(D)« 10" 
TTIME =TTIME* 107 [- 


T+ 
Go To UP eet te\ 
i410 9 


DOWN 


CLEAR LAL) tor bits tn 





Box 2 











DO |+07 
I= 1,MAA 
FO 1410 
L=1, LMA 





MASKSET 








| . Box 3 
TART okt DO 4 44q Com ube CA)= 
S L=\,LMAX BCL @ MASK Hd] 
A 1 
Shift CA) Left 4 = <i> —_ Switch 
Counts, NFACE =(A) a 
DOWN 
UP 
CALL 
ROUN DCL) 








Add one to 

NSIDE(NEACE). 
Add 8(C7,L) to 
ESTDECNFACE) 


tL: [L BCL)» MASK I3} 2 



















Te1bL 4 
TL Dye? 
Add one to NRBUL | _>O UP 
Add BG,L) to ERUL [<< 80, b) ~ THERMAL 
Set bit #2 of LB(L) 
DOWN 
Set bits of LBL) 









BC L)-ETHRE SH corresponding to bits 


| of Tu27 which ave set. 


Set bit #4 L_ 
of LBCL) | 


"On another page 





| | UP Go To 
Add B8C,L) tbo ETHERM , 


DOWN 
=. ae If 
£O ¢ =© Go To 
| DO 1426 Sule bk | bit #4 of Be 


Put L ‘in octal 
positions 6-4 
of LBCD). Put 







LBCT)*MASKI3}-IL27 
J in octal 
positions 6-4 
of LBCL) 





| 14.27 [A26 

Set bit +4 ey ZO. | Ald MASKSU | Go To 

of LB(T) boil to LBCL) 444 
[es 


we a Add MASKSA 
Set bit #4 | Go To 
of LBM) 44.9 * 
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1445] 


; Tio. 
Bethe @ of X(t) [LB(L)-MASKi] 


| If 
Go To #0 | Set bit #] 
of L&(1) 


=O 






Set MASKSU Tr. 2p 
bits of LB(LY [* Oo I Set MASKSA 
| bits of LBL) 









Do 1456 NVAC=O 
L=|,LMAX 


- toy 


Ac 
= 





Tf 
L *MASKSU a) 
Add one to 
NSUBST 


Box Add one to 
NYAC=NVAC+MAXLAT- LMAX (F@ 4-5G F 
INTERST 





| Box + 
. TEMP= Sum of ESIDE(I-6) | 
Sea ITEMP= Sumof NSIDECI-6) 


ECAL= ETHERM+ EBUL+ TEMP 
LMAXCAL= INTERST +NSAME + NSUBST+ NBUL + I TEMP 


PRINT: IX, ABCI), TY, ABC2), LZ, ABC), MAXLAT,A, 
SPHI, FMASS , ETHRESH, THERMAL | 


—=— =a wee 


L Ts LMAX-MAXLAT 


J 
[ 


YANGLE=. “32.aL PHA 
ZANGLE = +82. BETA 
















PRENT: IT, LULL, 
EPB, ALPHA, YANGLE, 
BETA, ZANGLE 





ITEMP=MAXLAT +1 | 






y 







PRINT: L, RET TT 
D ANU . wtS 
YEN TRY (L> MAXLAT),| | wyrs | a Faia 
} ~) 9 5 


ZENTRY(L- MAXLAT ) | 
\ for tion LMAX TSLICI-NAA), NTSCCI-NAA) 












PRINT: ECAL, 


TENERGY ,OENERSY, : 
ETHERM, CMAXC(D INTERST/2 


INTERST = 


PRINT: LMAXCAL , 
LMAX,NSAME ,NSUBST, 
TNTERST, NVAC 


PRINT: 
MISTAKE(I) 









PRINT: “ THE 
ENERGY CHECK 
OID NoT FAIL” 


2D PRINT: “NO INTERSTITIALS 
P| OR VACANCIES OCCURED ALONG 
AN EDGE OR CORNER” 


PRINT: NBUL, hee 
EBUL, NSIDE(i-<é) | 5 
ESIDE(|-6) 







PRINT: 
MISTAKE(2) 





Box 5 
PRINT: NSAME 
en LIVEM= 177777 70000777778 
M= |, LMAX 








Ti = 777 O000777I7T7II7F ITS 
J=KeL=) 
NUMBERS ()-10) = O 









Skip to next 
page o 
, output paper 


< 
_ 1474 = 


Skip one line on output 






NUMBERSC)=M 






PRINT: _ 
NUMBERS (TL) 
for IT=\, IT 
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ras PRINT: NVAC 


NUMBERSCIT)=L =O 
IT=IT+4I1 





Gt 


PRINT a 
NUMBERS (1-50) 


Box 
T 
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Box 6 


IT =| 


DO 1565 
L=1,MAXLAT 


IT=IT-| 


PRINT: 


NUMBERS(I-1T) 








Box / 
PRINT: NSUBST 


we DO: at 
= |, LMAX 


rL= [LBW MASKIS}F 
1 CALL INITIALCL) 
CALL. DEPTHPCL) 






PRIMM: L re, 
DEPTH, PERDIST 





DO \6§2l 
S17 tadas 








PRINT: 5, TEM, J=atof 


TEMP, Ic (I) 





TEM = A 0.5 
TEMP=J+0.5 


NUMBERS 1-30) = IC (I-30) DO 153I 


es) 









PRINT: LIT, 
eel (iseso) 





NUMBERS (1-20) = IC(31-50) CALL HIS To(20) Box 
he tes S 
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Box § 


PRINT: INTERST Ic(i- 100) = 0 


Go To #0 DO 1550 
iS 50 L=), LMAX 
> (A) = [LB(L)* MASKE9] 


TL= [LB(L)-MASKI3]+ 2” 7 


TE= DEPTH CALL INITIAL(L) 
TEM=PERDIST CALL DEPTHPCL) 


CALL INITIAL (LL) 








’ 






PRINT: L,LL,IL,TE,TEM, 





CALL DEPTHPCLL) DEPTH, PERDIST IS30 
PRINT: J, TEM, .. 2 DO 1556 IC C\- too)= 
TEMP, IC (1) naa T=1,30 Tcli-loof2 








NUMBERS (1-30) =1C (1-30) 
TL=-/ 


CALL DO 155% 

HISTO(30) T=], 20 : 
Box CALL N ot og PRET: 1, lt, € 

9 HISTO(I-20) wee TC (I+ 30) 
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2 


o> PRINT: NBUL asl 


Box 7 


=O 
Go To Bex JO 
Go To \57% 


DO 1569 ~ =6 
T= 1,30 IT=1T+I <y 
+O 









J=I-7 PRINT: L, B(-7,L), 
TEM=7-0.5 OCTAL(L), DEPTH, 
TEMPS I+0.5 PERPI'ST 


PRINT: J 
ee NUMBERS(1-30) =1C(1-30) 
rea J 8-7 





PRINT: rn 
1S 74 I,it,t, NUMBERS(I) = 
LC(L+30) _ _ICCT+30) 












PRINT: I, KETO), 
KET(Z+i), IC(I+50), 
for T=l,25 


Box | 
10 
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Skip to next 
page of out put 


CALL INITIAL(L) 
CALL DEPTHP(L) 
CALL ECOUNT(L) 


CALL 
HTSTO(30) 


DO 1574 
T=1,20 


NU MBERS(I-25) 


=IC(51-75) 





CALL 
H1STO(25) 





Rox 
DO 1599 <O_ ff PRINT: lO 
J —aleG NSIDEC\) 
>O 
PRINT: J Go To 
ZS jx : 
es NSIDECD) 1585 


NFACE =0 


DO 1595 
L= |, LMAX 


zo IC (\-100)=0 F 


NTS=0 










if 
Go To Go To =O | L) MAS A) 










NFACE = 


ee 


34 
iL B(L)- MASK 144-2 






NTS= 


[LB(L)- MASK25]> 


NFACE- J 





we 













PRINT: L, B(I-7,L), 
OCTAL(L), DEPTH, 
PERDIST, NTS 


CALL INITIAL(L) 
CALL CEPTHP(L) 
CALL ECOUNT(L) 











* Next page 
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a. a 





DO 1697 
I=|,20 


IT=1-l 






PRINT: L,TEM, 
TEMP, IC(1) 





[Pale = 7] 
TEMP O 0. 
= NUMBERS(I) = IC (I) 
CALL HISTO(30) | IL=-7 1696 


IL=O 


NUMBERS (I)= Re 1697 >| CALL HISTO(20) 


IC(I+30) 





iL, '@Gewee) 


0 
NUMBERS(Z)= IC(L+ 50) . a - 














PRINT: I, KET(I), 
KET(I+1), I1C(T +50) 
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Box I 
UP Go To 
120 


DOWN 


Go To 
avi 


Print out 


headings 


DO NNT 
L= 2. LMAX 












PRINT: L, 
B(\-7,L), OCTALCL) 










Print out 
headings 





<O 
20 
<0 


S kip one line 


PRINT: L, 117 on output 


Paper 


B(\-7,L), OCTALQ) 
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SUBROUTINE ROUND(L) Box I2 


IB(\-3)= TTEMP= 
B(I-3)/A +0.5 TBI) + 1B(2) + IB(3) 
IL= IB(3)* NPLANE+ EVEN 
RETURN 18(2)- NLINE + IBU)/2 +1 <> 
ODD 


NUMBERS(1) = 1B(3)* NPLANE + 1B(2)°NLINE +[TBC1) +1) 72 +1 
NUMBERS(2) = IB(3)- NPLANE +12B(2)-NLINE +[red\I)-)/2 +1 
NUMBERS(3) = 1B(3)* NPLANE + [1B (2)+ i]s NLINE + IBQ\)/2 + | 
NUMBERS(4) = I8(3) N PLANE +[I8(2)-\J- NLINE + r8(1)/2 +1 
NUMBERS (5)= [2 BC3)+ 1] NPLANE +2R(2)*NLINE+1B0)/2+1 
NUMBERS(6) =[TB(3)- Fs NPLANE+1B(2) NLINE+IRD/24+| 








TEMP= 10 bo 1340 =1=1,6 


DO) \as54,|_. CALL INITIAL(I1T) <O 
LTEMP= 1,3 TEM =0 


O 
2 Go To 
Add [xyz (x TEMP)- B(TTEMP,L) to TEM 


<O Tie oT 
arra> TEMP= TEM — 


= 
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= 


i 


~~ 


“ee 





SUBROUTINE INITIAL() eta 
<O TB(3)=(L- D/NPLANE 
1.6(2)= CQ)/NLINE 
TBCIYd= CA)°-2 
=1,3 
TBRA)=IBU)+1 


RETURN 










YYZ C\y= 0.0 
XY 2(2)= YENTRY(L-MAXLAT) 
KY 2(3) =ZENTRY(L-MAXLAT ) 
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( 


‘ 
€ 
° 
. 





SUBROUTINE ECOUNT(L) Boxes \4 and 15 


(sv TEMP = BC7,L)~KET(25) 
= a DO 20 <O 
seas, bb 


Ty. 50) + | 
es 


Box 
1S 






SUBROUTINE DEPTHPCL) 


<0 


PERDIST = 
SQRTF 4 


[e(2,)-xYzQ)]° 
+[e(,L)-xYz@} 





If 
PERDIST- 
As TEMP 


<O 


1640 
TC(r+30)= IC(T+30) + | (= 
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SUBROUTINE HISTO(NGROUPS) 


NTOTAL=0 

NUMBIG=0 
TOTALN = NTOTAL PRINT: 
BIGNUM = NUMBIG \. NTOTAL 


DO 113 iit 
J=1,NGROUPS 


Box |6 


Add NUMBERS(I) to NTOTAL 





. 


NUMBIG 


NUMBER= 100° NUMBERS(3) ZBIGNUM + 0.5 
TEMP = \100*NUMBERSCS)/ TOTALN 


ITEMP= NUMBER +4 vn 


DO N12 : 
ett. THISTO(T)= 3HIHX 
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= NUMBERS(CTI) 


DO III 
T=5, 104 


THISTO(Z)=3HIH 


Hi IT=JS4+1Ib 


ne PRINT: IT, TEMP 
with IHISTO as format 


=~ 


aes 


¢ P a 


i , : 
‘THE 


—_ - 





PROGRAM RON 









Write —\ 
3 times on 


Tape G 


Up 






DO 5 Skip to next page of output 
paper. Skip Cour li\MIeS on 
output paper 








Read one line Crecord) of 
Comment from Tape 2 and 
write it on output tape 










Read (Tape G): IX ,IN,7T2, ABCI-3), 
MAXLAT, A, SPHI, FMASS, ETHRESH, THERMAL 
MNP, {[B(2,5) T=1,3) J=1,MAXLATZ, NLINE,NPLAN 


Down 













Write on cutput tape: 
Lattice dimensions, 
number of atoms in 


lattice, lattice spacing 






LMAX= 
MAXLAT + MNP 


Okip four \tnes on 
out put paper 


Next 
Page 





constant, Sphere of 
influence, atomic mass , 









the two energy 
thresholds. 







2 


3 


wt 





Read T from — DX=T*A Set 
input tepe 2 RA=1/A N= | ve 
() ” <> Read IT 
=O 


and P(I-3) 
IN = upper half of IT DO 15 T= ABSEL 
IL = upper half of IT T=1,3 P(I)-B(I, 1L)] 


from Tape 6 
BC, IL) = PG) 
> 
X(I,N)= PCT) . 0 
NTOCN)= IT : 
<O 
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= 


‘ 1 ‘ 
’ e 
= . 
ee 
j o 
, é 





2 


M=0 DO 45 
00 45 ie : 
L=|,LMAX aA K=l,TTEMP 
LU= upper half =O II= lower half 
of NT(K) at Winck) 


#0 


Skip one 


line on 





*Preceeding page 
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APPENDIX V 
DESCRIPTION OF THE FLOW CHARTS 

The explanation given below is to be used with APPENDIX III 

(Program Listings) and APPENDIX IV (Program Flow Charts). 
PROGRAM MASTER 

BOX 1 

This box writes the title statement, "THIS IS PROGRAM MASTER" 
and as many comment cards (80 spaces per card) as desired. The DO 
loop is written for five comment cards at present, but can be readily 
changed if desired. 
BOX 2 

Reads input parameters for the run. Formats are given in the 
listing in APPENDIX III. 
BOX 3 

Changes lattice dimensions (IX, IY, IZ) into floating point format, 
and computes the physical dimensions of the lattice. These are then 
permanently found in AB (1-3). It must be kept in mind that IX must 
be ODD. 
BOX 4 

Computes NLINE, NPLANE, MAXLAT, LMAX (the initial LMAX) from 
the formulas given in APPENDIX II. 
BOX 5 

Computes initial lattice positions by use of subroutine SITE. 


Stores the site number and the fact that the site is occupied in the LB 
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ae 
eS tay" 
[| Ghee te Gag « « & cae 
tame | = at 8 a 
ee (le 





~ 
CO ASS 6 TST ee ET ay Coe oe oll 
Se - eee @ Oe — i 
Gut <8 oe on ee eee ee Oe) Fe cance 


ee ee 
- 








_——— 6 Oe ee Zee cee 
OT ona 
aT 

im pm ANNO HE :- |) rk game | 
ee SS ee ee eee oe ee ae 
A em: i ne Ree & mm 
iS 0 
‘- 

— = “ei 10 Se el) eee 
(ell ce ay ee Gemma ee 


- 








eS Ser CTee” ae eee ee 


= - | we art ea cee! @ oe 


array. Stores a -10.0 in each location of the TLAST array. This is 
necessary for collisions occurring in the first time step, otherwise 
these collisions would be ignored. 
BOX 6 

Writes initial lattice positions on Tape 6 in binary for eventual 
use by program RON. Also sets initial or zero values for the indicators 
shown. 
BOXES 7, 8, and 9 

Sets initial constant values, sets values of masks, and reads in 
the input parameters for the appropriate Born=-Mayer potential, either 
potential 1, 2, or 3. 
BOX 10 

This box includes the regeneration process. Jump Switches 1 & 2 
are raised, leaving Jump Switch 3 down. At the start of the next time 
step, the regeneration information will be written on Tape 7 in binary, 
and the tape will be rewound. To regenerate the run the program must 
be restarted, with all three jump switches in the raised position. After 
the program has started to read Tape 7 the jump switches may be reset. 
When the computer has finished reading Tape 7, the tape will be re- 
wound. Further regeneration tapes may be made after this point has 
been reached. The main DO loop over the time steps is also started in 
this box. Immediately after the start of a time step the current values 
of TENERGY, N, IENERGY, and ENERGY are written on Tapes 3 & 4. 


The output on Tape 4 is usually transferred to the typewriter to provide 
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ele 1 eked) 6 ace Sor ea eR Ce Tema al 
eS) A (OO Ts «(i ls a 


oo =: 






a ne 
nn: aaealetia 
1 A OM at eS ll 7 eT Pe 
og io li 

7 te 


wld oe ll 

[Te tio—— en” “Eee Cee a 
—, 
a Soe ee) © emp ee 0 cee © cel gee hae Aa, 
rs ce | 
2 ee ee i 
° ———~ & Se ow 

ees «/ .°- ee 














4 i 


| — 
4 mr et T Gece) ee, le eee 





— ee 
mere cleme,  aaeem eee: eer 


a) 
ee ee ee | 
“> 4 





‘= 
> ' ——_ -_ = © 0 ccc ‘ ~ 
Y cit ai «4 SS > 
ee, ‘ bo 7 
«? aoe 


* velo ete ot! Gael @ tale Sai 


a visual check on the validity of the run as it progresses. 
BOX 11 

All decisions concerning continuation or abandonment of the run 
are made in this box. A jump can be made to the output section, the 
normal output dump made, and a return to this box, or abandonment, 
accomplished by the outside control of the jump switches. The run 
will be abandoned, after an output dump, if the most energetic par- 
ticle has an energy less than THERMAL. However, the run will not be 
abandoned under this condition if other particles are to be shot into the 
lattice. 

Additional particles will be shot into the lattice, if the proper jump 
switch is set, at the start of the next time step that is one less thana 
multiple of LULL. For example, if LULL is 3, and the proper switch is 
set in time step 6 (or 7), an additional particle will enter at the start 
of time step 8. The switch must be left in the set position until after 
the particle has been read into the program. 

The input co-ordinates for the particle are given as - 1.0 for xX, 
and the actual point on the front face of the lattice we desire to hit. 
The angles ALPHA & BETA are given, as is the energy. This box then 
calculates the proper co-ordinates and velocity components to achieve 
the desired impact point. Since the new particle is the most energetic 
in the lattice, the Time Step Length is also recalculated at this time. 
BOX 12 


Here, the TSL is recalculated if the number of the time step is an 
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even multiple of JPB, and the first 30 values are placed in the TSLI 
table for output reference. 
BOX 13 
This is the only energy check made during the run, done every time 
step. The total energy (obtained as a sum of the individual energies) 
is compared with the original energy (the total amount shot into the 
lattice), and if they differ by more than 1%, the MISTAKE (1) error 
indicator is increased by one. 
BOX 14 
The value of ENERGY is recorded in the EMAX table for the first 30 
time steps, the total energy, TENERGY, is recalculated, and the LB 
array is reset and now shows that none of the atoms have been through 
the current time step. 
BOX 15 
The major DO loop on the particles is started here, with the LB 
array used to check certain conditions on each particle. If the atom: 
1) Does not have an energy greater than THERMAL, or 
2) Has left the lattice, or 
3) Has completed this time step, or 
4) Has been through INVAC and has not had a collision since, 
a jump is made to the end of the DO loop, and the next atom is con- 
sidered. 
Therefore, only those atoms are considered that satisfy all the 


following conditions: 
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1) The atom is a bullet, i.e. has more energy than THERMAL, 

2) The atom is in the lattice, 

3) The atom has not completed this time step, and 

4) The atom has not been through INVAC since its last collision. 
If all the preceeding conditions are satisfied, the initial values of the 
counters and lists shown in BOX 15, APPENDIX III are set to zero or other 
appropriate values. M, the fifth Index Register, is then set equal to L, 
the fourth Index Register. L, as such, is not used for the rest of the 
DO loop, until this pass through the DO loop is completed and a new 
L is chosen. 
BOX 16 

All atoms within a mathematical box with sides 3.02 (A) (centered 
at M) are placed on the NUM list here. Atoms are excluded if they have 
left the lattice or have been through the time step. If the list contains 
only one atom (M is placed first on the list), then a jump is made to BOX 
25 where another value of M is chosen from the LAST list. With other 
atoms in the box, indicators IN3 and NFACE are set to zero. 
BOX 17 
For every atom on the NUM list, the following parameters are com- 

puted with respect to atom M: 

1) The undeviated distance of closest approach between centers, 

called DSTANCE ( ), 
2) The time of closest approach between centers, also un- 


deviated, called TMIN ( ), 
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3) The time when the two atoms are separated by a distance of 
SPHI, called T ( ). This is the actual time a collision 
would occur. 

The distance is given in Angstroms, with the time in jiffys, relative 
to the start of the current time step. 

Various conditions are now imposed on these quantities to insure the 
validity of these presumed collisions. An atom must pass all the con- 
ditions to’insure that a collision will take place. The conditions are: 

=?) a x oe. 
2) DSTANCE ( ) § SPHI, 
3) T ( ) lies between + TSL of the start of this time step, 
4)T( ) is less than 0.99999 TSL. This eliminates collisions 
that end too far into the next time step. 
5) The TLAST for both atoms must be less than the absolute 
time of this collision. 

If all of the above conditions are satisfied, IN3 is increased by 
one. If one or more of the conditions is not satisfied, T ( ) is set 
equal to a large positive number (100 TSL serves this purpose in the 
program). 

It should be noted that DSTANCE (4), TMIN (4), and T (4) are the 
parameters associated with the atom that is number 4 on the NUM list, 
i.e. NUM (4), and so on. 


BOX 18 


The minimum time on the T ( ) list is found. All those on the NUM 
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list whose corresponding T's are equal to the minimum time, or within 
(CUTOFF) (TSL) of the minimum time are than placed on the KHIT list. 
Atom M is placed first on the KHIT list. At this point, if INPUT is zero, 
the KHIT list is duplicated on the LAST list. This occurs only when M 
equals L. INPUT is increased and the next box is entered if there is 
more than one atom on the KHIT list, i.e. another atom besides M. If 
there are no other atoms on the KHIT list, a jump is made to BOX 25, 
and M is set equal to the next atom on the LAST list. 

BOX 19 

This box and the next are concerned with the two body interaction 
itself, and the subsequent scaling methods for energy and velocity. 
These scaling methods are used whenever M collides with more than 
one atom. The pre-interaction velocities and energies are stored in 
the SAVE array. The COLD list is set to zero, and J is set to KHIT (JA). 
JA is the index used for the KHIT list in this box, M then hits the atoms 
in order of their placement on the KHIT list. JA is initially two, and 
increases up to and including NMAX. 

With J selected, LB (J) is checked, and if J is a member of an inter- 
stitial pair, NPAIR is set to 3. The indicator NPAIR is used in BOX 24 
to insure that the appropriate action is taken with respect to J's inter- 
stitial partner after the collision is over. 

Since the atoms on the NUM and KHIT lists do not have the same 
placement, a search is now made to find J on the NUM list. When this 


is done, we have the collision time, T, and the collision parameter, PAR. 
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The two atoms, M & J, are then moved ahead in time totime T. If M 
has already had one or more collisions, it must be moved back in time 
to the start of this time step, and then ahead to the start of this collision. 

The interaction subroutine is then called. Outputs are contained in 
the DEV array as the new positions, velocities, and energies of M & J. 

TLAST is now calculated for J, and J's co-ordinates and velocities 
are set to those of the DEV array. Notice that the TLAST for J does not 
include any o@leion of the time necessary for the collision. If a portion 
of this time were included, significant collisions would be eliminated. 
The COLD array (1-6) contains the sum of the final position and velocity 
co-ordinates for M as given by the subroutine for all of M's collisions. 
The total energy lost by M (ELOST) then equals the energy gained by 
this J plus that gained by the other J’s on the KHIT list. NHIT shows 
the actual number of atoms that have collided with M. 
BOX 20 

In this box, the scaling, by various means,of the velocities and 
energies of the atoms on the KHIT list is accomplished. 

The temporary quantity T2 is set equal to the previous energy of M 
minus the amount lost by M in all its collisions. 

If T2 is zero or negative’, the energy and velocities of M are set 
to zero. With T2 negative, the energies and velocities of all others 
on the KHIT list must be scaled to conserve energy. All atoms are 
scaled in the same proportion, the absolute direction of their velocity 


vectors is preserved, their energies being scaled by the square of 
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the scale factor used on the velocity components. 

In all cases, the spatial co-ordinates of M after the collision 
must be found. These are set equal to the numerical average of M's 
positions as computed by the subroutine. This is accomplished by 
simply dividing COLD (1-3) by NHIT. 

For a positive T2, there are two cases, the energy of M before the 
collisions above ETHRESH and below ETHRESH. 

When the original EY ETHRESH and there was only one collision, 
the co-ordinates of M are set to the DEV results, and a jump is made 
to the next box. If there was more than one collision, En is set equal 
to T2, and the velocities are scaled so that the components will be in 
the same proportion as the sum of the individual resultant velocities 
as computed by the subroutine. This essentially gives M a resultant 
direction equal to the vector sum of the directions M would have as a 
result of its collisions with the rest of the KHIT list. 

When the original E,, € ETHRESH, we assume that M has been 
oscillating about a site, and its site number is obtained from LB (M). 
Epis then set equal to T2, and the velocities are scaled so that the 
resultant velocity vector will point toward M’s oscillatory center. The 
scaling factor is (v2/ (AX2 +Ay? + Az2)) 1/2 This is multiplied by 
AX, AY, and AZ to obtain the correct components. This scaling auto- 
matically keeps We in proper proportion to Em 
BOX 21 


TLAST for M is calculated at the start of this box. M's TLAST 
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could not be calculated before this point since M had not yet parti- 
cipated in all its possible collisions. 

This box then insures that all atoms with energies less than 
ETHRESH have site numbers stored in LB, that they are placed on the 
LATER list if their energies are less than THERMAL, or that they are left 
alone if in between the two thresholds. 

If an atom has ES ETHRESH, its site number must be removed from 
LB (M), and the bit signitifying an occupied site must be removed from 
LB (site), unless the atom was once part of an interstitial pair. If it 
was part of an interstitial, then the partner section of LB must be re- 
moved from both LB's. 

If an atom should occupy a site and doesn't, then a site must be 
found. This is accomplished by going part way through INVAC, INVAC 
is exited by use of the MOVE indicator. Ifa site is calculated for M, 
and it is occupied, the information pertinent to an interstitial pair is 
recorded in LB. If the site is occupied by an interstitial, M is placed 
on the site anyway, and the fact that a "triple interstitial” was formed 
is recorded by increasing MISTAKE (5). EM is inspected at this point, 
and MISTAKE (3) is increased if Ea is negative. The symbolic language 
of this box may appear complicated, but comprehension should be 
facilitated by the flow chart. 

BOX 22 
All other atoms on the KHIT list are now subjected to the same pro- 


cess M underwent in the previous box, except that here a jump is made 
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to Statement 600 to enter INVAC. To preserve the values of I, J, andM, 
and also make M equal J, the values of I, J, and M are stored in MAD, 
NAVY, and LOCATE, respectively. These two boxes might, at some 
future stage of development be combined into one box covering the entire 
KHIT list. However, when this section was developed, the treatment of 
M and the others on the KHIT list was sufficiently different so that this 
combination would have been prohibitively difficult. Most of this dif- 
ficulty has disappeared as the program developed. 
BOX 23 

All atoms on the KHIT list are now backed up, time-wise and space- 
wise, to the beginning of the current time step. Each atom must be 
backed up using its particular T ( ), since they were not all hit at 
exactly the same time. A search through the NUM list accomplishes 
this, giving the correct value on the T list for the computation. 
BOX 24 

If one of the J's on the KHIT list was a member of an interstitial 
pair, then NPAIR is three, and a jump is made to Statement 236. The 
relations between atom M (the incoming atom), atom J (the member of 
the pair struck by M), and the partner, atom LL (the other member of the 
interstitial pair), are now considered. 

If J has an energy between the two thresholds (vibrating about its 
site), the information in LB about it must be retained. If the partner, 
LL, has been hit, and has an energy greater than ETHRESH, the MASK] 


portion of both LB'’s must be erased. If both have an energy ¢ ETHRESH, 
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nothing is done to either LB. If, however, J has an energy greater than 
ETHRESH and atom LL has not been hit, we must erase from LB (J) the 
site and partner number, erase from LB (LL) LL's partner number (which 
is J), and also move atom LL back onto the site. It would be a definite 
error to leave the partner, LL, on its previous split position, so it seems 
more realistic to move it back onto the site. 

BOX 25 

The interactions, scaling, and correction of LB have now been com- 
pleted for this round of collisions. A jump is now made to Statement 200 
in BOX 16. If M has no more collisions this time step, a jump back to 
this point will be made. If M does have further collisions, they are cal- 
culated and carried out. 

Upon entering this box, if IN3 is zero, then no collisions were made 
during the last pass through, and this M has completed the time step. 
The next atom on the LAST list, the INMIN th one,is made atom M. 
INMAX is the limit of this list. INMIN is increased by one, and a jump 
is made to Statement 200. When the LAST list has been exhausted, i.e. 
each one has been through the process until it will have no more colli- 
sions this time step, an exit to the next box is made. 

Any atoms on the LAST list that have been put on the LATER list to 
go through INVAC are not sent back through the collision process since 
they have MASK6 set. 

The name of the atom under primary consideration was changed from 


L to M in BOX 15 to prevent the above process from destroying the value 





of L, since Lis the index on the main DO loop over the atoms. 
BOX 26 

In a DO loop over the LAST list, only those atoms not on the LATER 
list are selected, and then moved to the end of the time step. This fact 
is then recorded in LB. The atom's energy is then compared with 
ENERGY, the largest one becomes the new ENERGY, and INERGY is 
changed if needed. 

BOX 27 

If NZ = 1, the LATER list is empty, and a jump is made around this 
box. If not, IT] is set equal to (NZ - 1), as the upper limit of the DO 
loop to follow. INVAC is then called successively for all of the LATER 
list, writing on Tape 6 in binary the particle number and time step, 
both contained in IU, and the X, Y, and Z co-ordinates. 

BOX 28 

In this DO loop on J, only those atoms on the LAST list that were 
not on the LATER list (and have consequently been through INVAC) are 
considered. 

All three co-ordinates are compared to the lattice dimensions. If 
an atom is outside the physical dimensions of the lattice, a value of 
NFACE is calculated. NFACE and the number of the time step of exit 
are then recorded in LB (M), the site number and possible interstitial 
partner are erased from LB (M), and M's number is erased from any 
interstitial partner's LB. The memory of occupancy is erased from the 


site if M is not an interstitial member, and the next atom on the LAST 
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list is considered. 
This is the end of the major DO loop on L, the next Lis chosen, 
or an exit is made to the next box if the DO loop on Lis completed. 
BOX 29 
The same procedure is followed here as was foliowed in BOXES 26 
and 28 above, except that atoms: 
1) That have been through the time step, 
2) Having an energy less than THERMAL, or, 
3) That are outside the lattice 

are not considered. 

This is the end of the major DO loop on N, the time step is completed, 
and either a new time step is started, or an exit to the next box and the 
output section is made. 

BOX 30 

Writes selected lattice and program parameters, anda — lattice 
dump on Tape 8 in binary for eventual use by the SLAVE program. Various 
BCD output is then written on Tape 3, and if the DO loop on N (the time 
steps) was completed, a (-1) is written in binary on Tape 6 and the pro- 
gram ends. If, however, the DO loop on N was not completed, a jump is 
made back to BOX 1l. There, a jump is made back to this box (with the 
resultant (-1) on Tape 6, etc.) if Jump Switch 3 is not set. In the normal 
operational mode, the DO loop is not completed, and the program itself 
terminates computation by a jump from BOX 11 to the start of this box, 


the output is written, the Jumps back to BOX 11 and back are made, and 
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the program ends. 

A jump around the BCD section of the output may be made by Setting 
Jump Switches 1 & 3. To facilitate resetting of the jump switches, Stop 
Switch 1 may be set. If this is done, the program will stop after writing 
Tape 8, and again before jumping back to BOX ll. 

Tapes 3 & 4, then, are used for BCD output, and Tapes 6, 7, and 8 
are used for binary output over the course of a run (including regeneration) . 
Since Tape 4 output is normally switched to the typewriter, only four out- 
put tapes are needed in the course of arun. The three binary tapes 
(6,7,8) will be rewound at the end of the run, and may be used without 
further action. An END OF FILE mark must be placed on the end of Tape 3 
before printing. 

SUBROUTINE SITE (M) 
BOK 31 

This short subroutine calculates the fixed point co-ordinates of 
the lattice site given as an input parameter. The co-ordinates are found 
in IB (1-3) after exit from the subroutine. 

SUBROUTINE INVAC (MJ) 

This subroutine, given the number of an atom, either places the atom 
on a vacant lattice site near it, or forms an interstitial pair with one of 
its nearest neighbors. 

BOX 32 
The initial indicator IN2 is set to zero, and the input atom is 


changed to number MM. The site number occupied by MM, and the fact 
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that MM is on the LATER list are now cleared from LB (MM), and the 
fact that MM has gone through INVAC is set in LB (MM). 

Atom MM's co-ordinates are now inspected to see whether or not MM 
is within one unit of a face of the lattice. If so, the number of this face 
is stored in IFACE (1-3), and the indicator NIF is increased. The index 
of IFACE is 1, 2, or 3, depending on whether the face is in the X, Y, or 
Z direction, respectively. 

The co-ordinates of MM are rounded off, and the fixed point co- 
ordinates of the nearest site (or non-site) are found. 

BOX 33 

We must now determine which of the two general cases we have, site 
or non-site, and if IFACE (1-3) is non-zero, which of the 12 special 
cases we are to consider. If the sum of the three fixed point co-ordinates 
is even, then they designate a lattice site, andJM is settol. If the 
sum is odd, then a non-site is designated, and JM is set to zero. The 
sum of the Y and Z co-ordinates is then checked. If odd, JL= 0, and if 
even, JL=-1. If the indicator MOVE is non-zero, a RETURN to the main 
program is executed. 

The formula for IL (the site number) results in an X co-ordinate 
one less or one more than an actual site for the —e case. The sum 
of JL and JM determine which case. If the sum is -l, then one is added 
rors. 

At this point, IL is always either: 


1) The site of the rounded off position of MM, or 
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2) The site with X co-ordinate one larger than MM's rounded-off 
position. 

The possible special case is now determined. If IFACE (1-3) = Q, then 
JB is set equal to the face number with JD = 1. If more than one IFACE 
is non-zero, MISTAKE (2) is increased, since this means a position along 
an edge or corner of the lattice. 

In either case, if IFACE (1-3) 3E0, then JA is non-zero, and JC is 
set equal to 1. If IFACE (1-3) is zero, then JC = JM (either 1 or 0). 
JB is set equal to the first non-zero IFACE or to zero if IFACE (1-3) 
is zero. JD is then set to 1 if a non-zero IFACE exists, and to zero 
otherwise. 

The case number (1-14) is now established by the formula: 
JW = §M*6*JD) + 1+ JB+JC. Cases 1 & 2 are the general cases. 
Cases 3-8 are special forms of case 1, with cases 9-14 as special 
forms of case 2. These special cases simply eliminate some of the 
nearest neighbors found in the general case from consideration. 

The lists MOON and NOON are set to zero, and the next box is 
entered. 
BOX 34 

The first section (down to Statement 2451) is a section of computed 
GO TO statements to insure the correct choice of the nearest neighbor 
formulas. The second section places the site numbers of the nearest 
neighbor positions on the MOON list, provided their numbers lie between 


one and MAXLAT. I is the index used for the MOON list, and is always 
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one greater than the number of sites on the list. 
BOX 35 

All of the nearest neighbor numbers found may not have been placed 
on the MOON list (possible negative numbers, or numbers greater than 
MAXLAT). JPMAX, therefore, is reset to (I-1). If the indicator MOVE is 
non-zero, a RETURN to the main program is executed. 

The LB for each atom whose site is on the MOON list is now in- 
spected. If the LB shows that the site is not occupied (see APPENDIX I), 
the site number is placed on the NOON list. JP is the index used for the 
NOON list. 

The MOON list, then, is a list of the sites of the nearest neighbors, 
and the NOON list is a list of those nearest neighbor sites that are vacant 
sites, with JK the number of vacancies on the NOON list. 

BOX 36 

If JK = 1, then only one vacancy exists, and atom MM is placed 
there, the site number, IL, is set in LB (MM), and the fact that the site 
is now occupied is stored in LB (IL). If IN2 is zero, a RETURN is then 
executed, otherwise a jump to BOX 39 (Statement 2800) is executed. 

If JK is zero, there are no vacancies and a prospective interstitial 
partner must be found. If IN2 is zero, a jump to BOX 38 is made to 
accomplish this. If IN2 is non-zero, then MM is the prospective inter- 
stitial partner of the original atom MJ, and a jump is made to BOX 39 to 


form the interstitial. 
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If JK > 1, then more than one vacancy exists. A jump is made to 
BOX 37, where one of the vacancies is chose, a retum to this box is then 
made, and the procedure given for JK = 1 is followed. 

BOX 37 

The (distance) @ between atom MM and the vacant sites on the NOON 
list is computed and stored in the DSTANCE array. The minimum Gee 
is then found, and the DSTANCE list is inspected to find the number of 
sites that are at this minimum distance from atom MM. lf there is only 
one, IL is set equal to that site, and a jump is made to Statement 2492 in 
BOX 36. If there is more than one, MM is placed on its own site, if its 
own site is one of those on the list. If not, then MM is placed on the 
first or second site on the NOON list, depending on whether the number 
of the time step is odd or even, respectively. 

The number of the time step is used here as a convenient random 
number generator, since the decision is made on oddness or evenness, 
rather than the actual value of the number itself. 

BOX 38 

This box is entered only from BOX 36, and then on the condition that 
IN2 is zero and that MM (MJ) has found no vacancies. A prospective 
interstitial partner must be chosen from the MOON list. 

The time of atom MM's (MJ's) closest approach to each of the sites 
listed on the MOON list is found and stored in the TMIN list. Any sites 


Occupied by interstitial pairs are not considered. 
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The minimum time of closest approach is then found, TMIN (KMIN), 

which locates the prospective interstitial site. MAX is derived from 

AX, AY, or AZ, whichever is greatest (Z\ X = X co-ordinate of MM 
minus X co-ordinate of the site chosen), and this becomes the prospective 
axis of the interstitial pair. The atom occupying this site is then found. 
Since it is possible for this atom to be pushed into any vacancy near it, 
and MM to fall into the vacancy thus created, this atom's nearest neigh- 
bors must be inspected. This is done by setting MM equal to the number 
of this atom, setting IN2 to 1, and jumping to Statement 2361 in BOX 32. 
BOX 39 

If the new MM (the original atom is still designated by MJ) does find 
a vacancy, it is placed there in BOX 36, a jump is made to this box, MM's 
now vacant site is named as the prospective site for MJ, a return to 
BOX 36 is made, and MJ is placed on its prospective site. A RETURN to 
the main program is then executed in BOX 36. 

If, however, the new MM did not find a vacant site, then a jump is 
made to this box (Statement 2720), and the interstitial is actually formed. 
The site of the interstitial is entered in each LB, together with the other 
atom's number. The fact that the site is occupied is then entered in the 
LB of the site, and a RETURN is executed. 

If the new MM equals MJ, then the atom has formed an interstitial 
with itself, and the error counter MISTAKE (4) is increased. The inter- 
stitial is formed so that MJ is on the same side as it was before forma- 


tion of the interstitial pair. 
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SUBROUTINE INTER (MQ, II, S) 

This subroutine, given the atom under consideration (MQ), the 
target atom (II), and the impact parameter (S, in meters) calculates the 
positions of the atoms at the end of the collision, their velocities, and 
energies and stores them in the DEV array. 

BOX 40 

Transforms atom MQ's position co-ordinates, velocity components, 
and energy and the angle between the radius vector from atom II to 
atom MQ and MQ's velocity into the target (atom II) frame of reference. 
Sets CPA = SPHI as the initial value for the CPA iteration. 

BOX 41 

Computes CPA by an iterative process. The iteration is terminated 
when two succesSive values differ by less than 1074 %, or when any 
value is greater than the preceeding value. 

BOX 42 

Computes the time required for the interaction and the total angle 
of deviation. The integrals are solved by numerical integration using a 
Gauss-Legendre four point quadrature. 

BOX 43 

Computes the composite factors and trigonometric functions that are 
used more than once in the transformation to the laboratory reference 
frame. 

BOX 44 


Computes direction cosines, relative to the laboratory frame of 
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reference, of a line in the plane of the interaction that is perpendicular 
to the line between the two atoms at the beginning of the interaction. 
BOX 45 

Computes the position co-ordinates, velocity components, and 
energies of both atoms and stores them in the DEV array. See APPENDIX I, 
DEV for the actual location. 

PROGRAM SLAVE 
Bem | 

Reads binary input from Tape 8. "SET SLJ SWITCH 2 UP IF NEEDED" 
is written on Tape 3 (normally transferred to the typewriter). As many 
comment cards (80 spaces per card) as desired (6 at present) are read 
from Tape 2 and written on the output tape. Reads input parameters for 
the run from Tape 2. Writes on the output tape "SELECTIVE JUMP SWITCH 
NUMBER TWO WAS UP" if Jump Switch 2 is used for the run. 

IHISTO (1-4) and IHISTO (105) are defined, the IHISTO array then 
serves as the format specification for the histographs. The constants 
IZ2E15, I12E24, F2E24, and 12E30 are set, and the masks for use with the 
LB array are defined. 

BOX 2 

If NTS (called N in PROGRAM MASTER, the number of time steps 
completed) is zero, then NTS is set equal to MNTS; other constants and 
indicators are then set. The times in the TSLI array and the time TSL 
are converted from jiffys to seconds. Various portions of the LB array 


are then cleared, depending on the settings of the jump switches. 
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BOX 3 
In this box all atoms are placed in one of the following classifica- 
tions: 
1) Atoms occupying their original sites, 
2) Atoms occupying sites other than their original sites, 
3) Atoms that are parts of interstitial pairs, 
4) Atoms whose energy is greater than ETHRESH, and 
5) Atoms outside the lattice. 
Program MASTER has also calculated the face of exit for those atoms that 
have left the lattice. 
The number oz atoms and the total energy of those atoms is calculated 
for each classification, as is the number of vacant lattice sites. 
BOX 4 
The total number of atoms and their total energy is calculated. 
This serves as an additional total energy check at the end of the compu- 
tational run. These numbers and various other data (see APPENDIX III) 
are written on the output tape. The number of atoms and their total 
energy is written on the output tape for each classification listed in 
BOX 3. 
BOX 5. 
The number of atoms occupying their original sites is written on the 
output tape. An array of these atoms is then written, in order, on the out- 
put tape. Atoms not on their original sites are represented by zeros in 


this array. 
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BOX 6 


The total number of vacant lattice sites and the numbers of these 
sites are written on the output tape. 
BOX 7 

The total number of replacements is written on the output tape. 

The number of the atom, the number of its final site, its penetration, 

and its radial penetration are written on the output tape for each of 

the replacement atoms. Histographs are then compiled and written on the 
output tape for the penetration and the radial penetration distributions 

of the replacement atoms. 

BOX 8 

The total number of interstitial pairs is written on the output tape. 
The number of each atom in a pair, its partner, the site they occupy, the 
penetration of each atom in the pair, and the radial penetration of each 
atom in the pair is written on the output tape. This is written twice for 
each interstitial pair since there are two atoms in the pair. Histographs 
are compiled and written on the output tape for the penetration and the 
radial penetration distributions of all interstitial pair members. 

If Jump Switch 2 was set for the run, the interstitials are determined 
by inspecting the LB array given by program MASTER. If, however, Jump 
Switch 2 was not set, the program disregards the LB array and calculates 
the interstitial pairs in the lattice. 

BOX 9 


The total number of atoms with energy greater than ETHRESH 
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(knock-ons) is written on the output tape. The number of each atom, its 
co-ordinates, velocity components, energy, its LB (in octal), its pene- 
tration, and its radial penetration are written on the output tape for each 
knock-on atom. Histographs are then compiled and written on the out- 
put tape for the penetration, radial penetration, and energy distributions 
of the knock-on atoms. 
BOX 10 

The total number of atoms that exited through a lattice face is 
written on the output tape. The number of the atom, its co-ordinates, 
its velocity components, energy, LB (in octal), its penetration, radial 
penetration, and time step of exit are written on the output tape for 
each atom exiting through a face. 

The entire box is then repeated for each of the six faces of the 
lattice. 

BOX 11 

If Jump Switch 1 is set, a jump is made to the end of the program, 
otherwise a general lattice dump is written on the output tape. The dump 
consists of the atom number, co-ordinates, velocity components, energy, 
and the LB contents (in octal) for each atom. 

SUBROUTINE ROUND (L) 

BOX 12 

This subroutine computes the site closest to atom L. The atom's 
co-ordinates are rounded off, and the fixed point co-ordinates of the 


atom are stored in the IB array if the co-ordinates represent a site. 
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If the fixed point co-ordinates do not represent a site, the six nearest 
neighbor sites are found. The site that is closest space-wise to the 
atom is chosen as the site for atom L. Once a site is chosen, a RETURN 
to the main program is executed. 
SUBROUTINE INITIAL (L) 
BOX 13 
If atom L is one of the original lattice atoms, this subroutine 
computes its initial fixed point co-ordinates and stores them in the IB 
array. If atom L was shot into the lattice, the point of entry is stored 
in the IB array. A RETURN to the main program is then executed. 
SUBROUTINE ECOUNT (L) 
BOX 14 
The energy of atom Lis inspected. This energy will fall between 
two of the values on the KET list, say the sixth and the seventh values 
on the list. Atom L's energy is then greater than the sixth (Ith) value 
on the list. The histograph counter IC (I + 50) is then increased by one. 
A RETURN to the main program is then executed. 
SUBROUTINE DEPTH (L) 
BOX 15 
The penetration (DEPTH) is computed as the final X co-ordinate of 
atom L minus the initial X co-ordinate. A value of I is determined so 
that DEPTH lies between A (I -7.5) and A (I -6.5). Histograph counter 


IC (I) is then increased by one. 
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The radial penetration (PERDIST) is then calculated from L's initial 
and final positions. A value of I is determined so that PERDIST lies be- 
tween A (I -1) and A (I). Histograph counter IC (I + 30) is then increased 
by one. A RETURN to the main program is then executed. 

SUBROUTINE HISTO (NGROUPS) 
BOX 16 

The largest number and the total of NUMBERS (1 -NGROUPS) are 
determined. The total is written on the output tape. The percentage of 
the total is written on the output tape for each entry on the NUMBERS 
list, followed by a number of X's corresponding to the percentage this 
entry is of the largest number in the table, i.e. the largest number will 
write out 100 X's, and a number half that large will write out 50 X's. 

A RETURN to the main program is then executed. 
PROGRAM RON 

This program writes on the output tape the successive positions of 
every atom in the lattice (or that was shot into the lattice), if the suc- 
cessive position of an atom is significantly different from the previous 
position. The number of the time step is also written when an atom's 
co-ordinates are written. After writing all the successive positions of 
an atom, the program commences writing the positions of the next atom 


that has moved significantly until all atoms have been considered. 
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